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INTRODUCTION 


Because  it  is  the  west  coast  satellite  launch  center  for  DOD, 
Vandenberg  AFB  supports  on-going  improvement  of  atmospheric 
dispersion  models  used  in  risk  assessment  for  potential  "hot"  and 
"cold  spills"  of  stored  rocket  propellents  and  oxidizers  or  those 
in  transport  or  transfer  around  the  active  launch  sites. 

Potential  "hot  spill"  assessments  are  now  treated  with  the  REED-M 
model.  "Cold  spills"  use  the  Ocean  Breeze/Dry  Gulch  model  (to  be 
replaced  by  the  AFTOX  Gaussian  plume  model) .  More  defensible  "hot" 
and  "cold  spill"  models  have  been  mandated  for  emergency  response, 
normal  launch,  and  routine  venting  purposes.  Thus,  together  with 
the  other  ranges,  Edwards  AFB  and  Cape  Canaveral,  Vandenberg  is 
actively  considering  a  number  of  diagnostic  (time  invariant)  and 
prognostic  (time  varying)  models  for  on-base  operational  use. 

The  creation  of  timely,  accurate  predictions  using  prognostic 
models  strains  the  limits  of  existing  computers.  So  the  Vandenberg 
risk  assessment  community  has  deemed  diagnostic  model  upgrading  the 
initial  priority.  As  available  computer  power  increases,  the  focus 
will  shift  to  selecting  and  tailoring  improved  prognostic  models. 

Within  this  two-tiered  approach,  step  one  is  to  evaluate  end  user 
needs  at  the  Vandenberg  launch  pads  and  storage  facilities.  This 
was  done  by  survey,  with  input  from  other  test  ranges.  Next  was 
preliminary  survey  of  all  available  existing  dispersion  models. 
Twelve  candidate  models  were  discussed  at  the  USAF  Toxic  Release 
Advisory  Group  semi-annual  meeting.  Further  assessment  narrowed 
the  group  to  the  following  serious  contenders:  NUATMOS/CITPUFF, 
CALMET/CALPUFF,  PGEMS,  WOCSS/MACHWIND/Adaptive  plume,  LINCOM/ 
RIMPUFF/HEAVYPUFF,  MATHEW/ ADPIC .  The  following  evaluates  these  six 
diagnostic  atmospheric  dispersion  modeling  groups  on  the  basis  of 
technical  merit  and  suitability  for  the  Vandenberg  environment. 
This  is  being  pursued  as  part  of  Phase  I  of  the  diagnostic  model 
evaluation,  testing,  validation,  and  installation  effort. 

We  attempt  here  to  recommend  an  optimal  modeling  grouping, 
preliminary  to  detailed  code  testing,  operational  simulation,  and 
the  validation  and  installation  phases. 

We  begin  with  lists  of  important  modeling  issues  and  models, 
specific  model  features,  Vandenberg  requirements,  and  figures 
assessing  the  current  cost  and  availability  of  dispersion  modeling 
resolution  and  computer  hardware.  We  then  assess  the  degree  of 
fulfillment  of  Vandenberg  requirements  and  modeling  issues  on  the 
basis  of  available  documentation.  The  format  consists  of  a 
narrative  description  and  evaluation  of  six  modeling  groups  and 
summary.  The  appended  reference  list  is  lengthy,  but  includes  much 
of  the  standard  literature  considered  useful  in  assessing  the 
current  state  of  atmospheric  dispersion  modeling. 
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MODEL  ISSUES 


1.  Cost  and  Ease  of  Procurement  and  Use,  Life-Cycle  Utility 

a.  Site  license  availability 

b.  Right  to  modify 

c.  Technical  support 

d.  Right  to  upgrades 

e.  Total  hardware/software  cost  including  transition  costs 
from  current  system 

f.  Estimated  time  and  effort  for  validation 

g.  Time  to  operational  status 

i.  Current  model  readiness 

ii.  Effort/benefit  analysis 

h.  Estimated  training  time  and  user  salary 

i.  Estimated  time  to  obsolescence 

2.  Military,  Regulatory,  and  Legal  Requirements 

a.  Military  mandates 

b.  Beyond  Base  regional/state/federal  regulatory  requirements 

c.  Legal  def ensibility 

4.  Input  Data  Requirements,  adequacy  of  current  measurements 

a.  Towers 

b.  Sodar 

c.  Rawinsonde 

d.  Radar 

e .  Buoy 

f.  Radiation/Cloud  cover 

g.  Soil  parameters 

h.  Variable/homogeneous  initialization 

5.  Domain  issues 

a.  Domain  size 

b.  Grid  spacing 

c.  Domain  type  constraints:  hills,  land/water,  buildings 

d.  Surface,  Lateral,  and  Top  Boundary  conditions 

6.  Windflow  Model 

a.  Initialization  procedures  and  objective  analysis 

b.  Stable/Neutral/Unstable  Physics 

i.  Turbulence  closure  order 

ii.  Moisture  treatment 

iii.  Boundary  layer  parametrizations 

c.  Grid  solver,  e.g.  finite  difference/element/spectral 
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Diffusion  Model 


a.  Source  parameters 

b.  Chemistry 

c.  Type  of  physics,  e.g.,  Gaussian  plume,  puff,  particle 

d.  Stable/Neutral/Unstable  Physics 

e.  Dense/Neutral/Buoyant  gas  treatment 

f.  Numerical  solution  procedures 

Fire  and  Explosion  Model 
Risk  Assessment  Model 

Total  system  issues 

a.  Comprehensiveness 

b.  Internal  consistency,  module  sophistication  balance 

c.  Modularity/Upgradability 

d.  Portability 


MAJOR  MODELS  EVALUATED 


1 .  CALMET/CALPUFP 

Sigma  Research  (J.S.  Scire,  D.G.  Strimaitis,  R.J.  Yamartino) 
for  California  Air  Resources  Board  regional  air  quality 


2.  POEMS 

Battelle  Northwest  Laboratory  (K.J.  Allwirie,  C.D.  Whiteman,  V. 
Ramsdell) 

for  Pacific  Gas  and  Electric  Co.  at  Diablo  Canyon  Nuclear 
Power  Plant 

3.  NUATMOS/CITPUFF 

USDA  (D.G.  Ross,  D.G.  Fox) 

for  U.S.  Environmental  Protection  Agency 

4.  WOCSS/Adaptive  Plume 

SRI  Int.  (F.L.  Ludwig,  Roy  Endlich) 
for  SRI  Int. 

5.  MACHWIND 

U.S.  Army  Atmospheric  Sciences  Laboratory  (R.  Meyers) 

G .  LINCOM/RIMPUFF/HEAVYPUFF 

RISO  National  Laboratory/Naval  Postgraduate  School 

(I.  Troen,  G.  Lai,  T.  Mikkelsen,  M.  Nielsen) 

for  USAF  Space  and  Missile  Systems  Center  Vandenberg  AFB 

THC  dispersion  model 

7.  MATHEW/ADPIC 

Lawrence  Livermore  Laboratory  (M.  Dickerson,  R.  Lange,  D.  Ermak) 
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VANDENBERG  AFB  SPECIFIC  REQUIREMENTS 


CALMET/ 

PGEMS 

WOCSS/ 

LINCOM/ 

NUATMOS 

MATHEW/ 

CALPUFF 

Adaptive/ 

RIMPUFF/ 

CITPUFF 

ADPIC 

MANDATORY  REQUIREMENTS: 

MACHWIND 

HEAVYPUF 

a . 

Treats  complex  terrain 

3 

3 

3 

3 

3 

3 

ta. 

Gives  surface  footprints 

3 

3 

3 

3 

3 

3 

c . 

Includes  all  met  data 

3 

3 

3 

3 

j 

3 

d. 

Handles  50x50km  domain 

3 

3 

3 

3 

? 

3 

e . 

Modeled  winds  near  measured 

3 

3 

3 

3 

7 

3 

f . 

Handles  gas/liq  puff/plumes 

0 

0 

0 

1 

0 

1 

g- 

"  surface/elevated  sources 

3 

3 

3 

3 

3 

3 

h. 

"  spills  by  v;eight,  volume. 

0 

0 

0 

0 

0 

0 

wetted  area,  &  flow  rate 

i . 

"  dense  gas  &  VBG  chemistry 

0 

0 

0 

1 

0 

2 

j  • 

"  on-line  data  &  gives  THCs 

3 

3 

3 

3 

0 

3 

k. 

"  multiple  THC3 

3 

3 

0 

3 

9 

3 

1. 

Outputs  THCs  w/i  5  min. 

0 

3 

3 

3 

9 

3 

m. 

Allows  user  overrides  in  GUI 

0 

1. 

5 

0 

3 

9 

3 

SIGNIFICANT  REQUIREMENTS: 

a , 

Adequate  grid  resolution 

0 

0 

2 

2 

9 

2 

b. 

Treats  variable  Zi 

2 

2 

2 

1 

9 

1 

c . 

"  all  PGFs  &  stagnation 

2 

2 

2 

1 

2 

1 

d. 

"  wet  and  dry  deposition 

2 

2 

0 

2 

? 

2 

e . 

Graphical  &  tabular  output 

2 

2 

2 

2 

9 

2 

f . 

Multi-level  output 

2 

2 

2 

2 

2 

2 

g- 

Treats  vector  shear 

2 

2 

2 

2 

9 

2 

h. 

Source  code,  docs  available 

2 

2 

2 

2 

? 

2 

i . 

Site  license  &  mod  rights 

2 

2 

2 

2 

9 

1 

DESIRED  REQUIREMENTS: 

a . 

Modular  code 

1 

1 

1 

1 

? 

1 

b. 

MARSS  compatible 

1 

1 

1 

1 

1 

1 

c . 

Flexible  chem  module 

1 

0. 

5 

0 

0 

0 

0 

d. 

Treats  coastline  &  cloud/ 
clear  interfaces 

0 

0 

0 

0 

9 

0 

e . 

Operates  w/i  current  MARSS 

0 

0 

0 

0 

0 

0 

f . 

building  wake  effects 

1 

1 

0 

0 

0 

0 

Total  score 

44 

48 

42 

50 

7 

50 

out  of  Oj  possible  points 


CAVEAT;  Detailed  decisions  compel  hard  thinking  which  cannot  be  distilled  into 
a  simple  score.  That  the  above  scores  are  mainly  congruent  with  our  basic 
impressions  is  largely  fortuitous.  From  the  above  comparison  (allowing  up  to 
3  points  for  mandatory  features,  2  for  significant,  and  1  for  desirable,  and 
fractional  scores  for  partial  fulfillment),  it  seems  that  LINCOM/RIMPUFF  and 
MATHEW/ADPIC  display  the  most  elements  desired  for  a  Vandenberg  atmospheric 
dispersion  model.  Actually,  since  windflow  and  diffusion  are  usually  separable, 
a  combination  of  LINCOM,  MACHWIND,  and  ADPIC,  or  either  RIMPUFF  or  CALPUFF  may 
be  the  optimal  overall  choice,  if  obtainable.  For  our  specific  reasons,  see  the 
following  text.  We  also  do  not  feel  that  PGEMS  comes  as  close  to  the  top  models 
as  is  indicated  by  its  score. 
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COMMENTS  ON  FIGURES 


The  following  figures  are  intended  to  illustrate  the  cost  of 
increasing  grid  resolution  for  dispersion  models  in  terms  of 
computer  hardware  as  projected  for  calendar  year,  1993.  Figure  1 
compares  the  cost  of  state-of-the-art  representative  PCs, 
workstations,  mainframes,  and  supercomputers  against  their 
computational  speed,  measured  in  units  of  (MFLOPS)  one  million 
floating  point  (real  number)  instructions  per  second.  Four  points 
may  be  made.  1)  The  sample  older  installed  systems  at  Vandenberg 
are  less  cost  effective  by  one  to  two  orders  of  magnitude  than 
state-of-the-art  systems.  2)  The  speed/cost  ratio  increases  by 
perhaps  an  order  of  magnitude  as  modern  systems  get  smaller,  due  to 
simplifications  in  job  control  scheduling,  peripherals  designed  for 
smaller  numbers  of  users,  and  efficiencies  of  production  scale.  3) 
Workstation  and  PC  systems  are  about  to  drop  below  $1,000  per 
MFLOP.  4)  The  ultimate  potential  of  massively  parallel  systems 
exceeds  other  approaches,  but  current  combinations  such  as  Intel's 
128  i860  CPUs  have  not  broken  away  yet  from  the  overall  curve,  nor 
are  they  yet  as  fast  on  real  world  problems  as  older  vector 
supercomputer  designs  such  as  the  CRAY  C-90.  This  may  be  due  to 
immaturity  in  hard/software  systems  and  compiler  design,  as  well  as 
inadequate  communications  bandwidths  between  processors. 

Figure  2  illustrates  the  effects  of  computer  speed  on  allowed  model 
resolution  for  sample  operational  diagnostic  and  prognostic 
dispersion  models,  as  well  as  a  sample  research  grade  large  eddy 
simulation  (LES)  windflow  model.  We  see  that  1)  Diagnostic  models 
require  less  computer  power  than  prognostic  ones,  while  LES  models 
are  very  computer  intensive.  2)  The  speed/resolution  slope  is 
steeper  for  prognostic  models  because  for  prognostic  models  greater 
spatial  resolution  also  requires  greater  time  resolution,  while 
diagnostic  models  disregard  time  variations  entirely. 

Figure  3  combines  figs.  1  and  2  to  suggest  the  hardware  cost  of 
resolution  for  the  same  sample  models.  Diagnostic  models  can  or  at 
least  should  run  effectively  on  systems  in  the  one  to  ten  megaflop 
speed  range.  This  means  that  efficient  operationally  oriented 
diagnostics  models  should  be  able  to  run  effectively  on 
workstations  priced  at  $10,000  to  $60,000  and  available  now;  while 
for  the  forseeable  future  LES  models  will  still  require  the 
resources  of  vector  or  massively  parallel  supercomputers  above  the 
$1,000,000  price  point  and  will  remain  research  tools.  This 
suggests  that  emergency  response  systems  no  longer  need  to  depend 
on  tenuous  linkages  with  remote  supercomputers.  A  second  point  can 
be  illustrated  in  more  detail  by  inspecting  the  computational 
configuration  ASTeR  has  suggested  for  operational  use  of  RAMS  at 
Cape  Canaveral  (Lyons,  et  al.,  1992) 
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GRID  NO. 

X 

Y 

2 

AX(km) 

TIME  STEP (sec) 

GRID  PTS 

1 

38 

X 

38 

X 

25 

60 

90 

36,100 

2 

34 

X 

38 

X 

25 

15 

30 

32,300 

3 

37 

X 

37 

X 

18 

3 

12 

TOTAL 

24 , 600 

93,000 

Since  the  spatial  nesting  ratios  are  already  quite  strained  within 
this  configuration,  running  such  a  model  at  500  meter  resolution 
over  Vandenberg  might  require  more  than  one  additional  level  of 
grid  nesting.  To  account  for  local  vertical  accelerations  in 
terrain,  at  least  the  inner  two  grids  must  operate  non- 
hydrostatically ,  perhaps  with  second  order  turbulence  closures, 
while  the  inner  grid  uses  a  2  -  3  second  time  step.  Such 
modifications  would  require  more  than  an  order  of  magnitude  more 
computer  power  than  is  now  available  on  workstations.  The  cost  of 
computing  power  has  dropped  by  about  a  factor  of  two  every  two 
years  for  the  past  decade.  So  it  seems  for  Vandenberg  that 
operational  fine  scale  prognostic  modeling  on  workstations  may  be 
some  years  away.  Massively  parallel  systems  may  eventually 
accelerate  computer  power  gains,  but  the  timing  does  not  seem  near 
term  and  beyond  that  is  difficult  to  gauge. 

In  the  interim  and  perhaps  beyond  it,  there  exists  a  niche  for  fine 
scale  diagnostic  models.  In  fact  a  hybrid  solution  consisting  of 
a  prognostic  model  which  supplies  coarse  resolution  mean  winds  to 
a  fine  scale  diagnostic  model  is  quite  defensible  and  could  be 
tuned  to  run  in  seconds.  Such  a  hybrid  could  conceivably  supply 
250  meter  resolution  winds  based  on  currently  available  hardware. 
As  available  computer  power  increases,  the  resolution  of  such  a 
hybrid  may  fall  to  the  100  meter  level  where  some  believe  the  point 
of  diminishing  returns  begins. 
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MEGAFLOPS 


COMPUTER  COST  vs  SPEED  (1993  $) 


CrayCaO 


10  "  10 
DOLLARS 


FIG.  1 
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MEGAFLOPS 


MODEL  RESOLUTION  VS  COMPUTER  SPEED 


No.  of  GRID  PTS 


FIG.  2 
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3,000,000 


(ote  that  the  immediately  following  CALMET/CALPUFF  review  based  on 
'.<=  user's  manuals  #A025~1,2  and  personal  communications  with  the 
-aiif.  Air  Resources  Board  and  University  of  Calif,  at  Davis  is 
lengthier  than  reviews  of  some  of  the  other  models.  If  the 
critique  seems  more  detailed,  it  is  often  because  the  comments  also 
apply  to  many  of  the  other  models  and  diagnostic  models  in  general. 

CALMET/CAI.PUFF  (J.S.  Scire,  et  al.): 


INTRO: 

CALMET/CALPUFF  combines  either  a  diagnostic  or  prognostic  wind  flow 
model  with  fairly  sophisticated  puff  dispersion.  The  output  is  now 
tabular,  but  UC  Davis  is  under  contract  with  the  California  Air 
Resources  Board  (CARB)  to  evaluate  CALMET/CALPUFF  and  apparently 
also  provide  an  output  format  using  NCAR  graphics.  A  full 
graphical  user  interface  is  not  planned.  A  DEC  VMS  version  of  the 
Fortran  code  exists  and  CALMET/CALPUFF  has  been  run  on  a  MicrcVax 
II,  as  well  as  DEC  3000/200  and  5000/200  Stations,  so  MARSS 
compatibility  is  probably  not  an  issue.  CALMET/CALPUFF  has  a 
limited  number  of  user  override  options. 

INPUT; 

In  its  current  form  CALMET  proceeds  from  domain  average  winds  which 
can  be  specified  from  a  sounding.  It  also  uses  surface  tower  data, 
and  a  single  domain-wide  vertically  averaged  temperature  lapse 
rate.  CALMET  evaluations  in  progress  at  UC  Davis  indicate  that  the 
model  needs  to  be  altered  to  accept  conditions  where  data  is 
missing  or  when  the  available  sounding  does  not  reach  the  specified 
top  of  the  domain  (default  1200  meters) . 

PHYSICS: 

CALMET  computes  step  1  and  step  2  wind  fields  with  step  1  having 
options,  a  or  b.  In  option  a,  the  time  interpolated  sounding  or 
user  specified  domain  average  wind  is  adjusted  empirically  for 
terrain  effects  due  to  kinematics,  slope  flow,  and  blocking.  The 
blocking  occurs  according  to  a  critical  streamline  height  (CSH) 
defined  by  a  local  Froude  number  criterion,  Fr  =  U/NAh  <  1,  (See 
WOCSS  review  for  more  detail) .  Here  U  is  local  wind  speed.  Ah  is 
the  difference  between  local  terrain  height  and  the  highest  height 
within  the  local  area.  The  Brunt-Vaisala  frequency  is  defined  as, 
N  =  [  (g/0)  (30/3z)  ] '''^,  where  g  is  gravitational  acceleration  and  0 
-  T(Po/p)’^^P  is  potential  temperature.  p  and  Po  are  pressure  and 
surface  pressure,  R  is  the  gas  constant  and  Cp  is  the  heat  capacity 
of  air.  The  significance  of  the  CSH  is  that  below  it,  flow  is 
thought  to  diverge  tangentially  around  the  hill,  while  above  it 
flow  proceeds  over  the  hill.  The  Froude  criterion  can  be  applied 
straightforwardly  to  an  isolated  hill  but  may  be  ambiguous  in 
general  complex  terrain.  The  code  should  be  inspected  carefully  to 
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assure  that  the  Froude  concept,  will  apply  in  a  reasonable  manner 
for  the  possible  wind  conditions  and  Vandenberg  terrain. 

In  practice  the  CSH  concept  is  widely  accepted  and  used  in  most 
diagnostic  models  which  operate  for  stable  conditions.  It  is 
consistent  with  the  results  of  laboratory  expex'iments,  using  a 
number  of  hill  shapes  (Hunt  et  al.  1988b;  Meroney,  1990) .  The  CSH 
is  based  on  the  assumption  of  a  kinetic  versus  potential  energy 
balance.  I.e.,  an  air  parcel  completely  converts  kinetic  to 
potential  energy  in  rising  to  some  maximum  attainable  height 
indicated  by  the  critical  streamline.  Though  the  values  calculated 
for  CSH  have  been  consistent  with  data,  the  underlying  Froude 
concept  is  itself  physically  suspect.  For  one  thing  it  implies  a 
zero  windspeed  at  the  CSH,  as  well  as  local  minimums  atop  hills; 
both  of  these  features  are  contrary  to  all  experimental  data.  The 
Froude  number  criterion  ignores  the  role  of  horizontal  pressure 
gradients  on  windspeeds.  By  considering  hydrostatic  and  Bernoulli 
constraints  on  flow  in  a  Lagrangian  frame.  Smith  (1990)  has  shown 
that  pressure  gradients,  rather  than  potential  energy,  control  wind 
speeds  and  these  gradients  vary  with  hill  steepness  and  aspect 
ratio.  So  the  success  of  the  Froude  criterion  may  be  fortuitous. 

Eventually,  we  would  prefer  to  see  models  for  stable  conditions 
which  are  consistent  with  Smith's  more  physically  realistic  basis. 
In  particular.  Smith  has  shown  through  linear  approximation  that 
there  will  be  two  regions  of  minimum  v/indspeed,  one  above  the  hill 
top  at  height,  z  =  4 . 0  U/N,  and  one  near  the  hill  surface  along  the 
upwind  slope.  The  height.  Ah,  of  the  near-surface  velocity  minimum 
varies  from  roughly  U/N  to  3U/N,  as  the  hill  aspect  ratio  varies 
from  flow  perpendicular  to  a  thin  ridge  to  flow  over  a  two- 
dimensional  hill.  For  a  radially  symmetric  hill  the  height 
estimate  is  Ah  =  1.3U/N.  So  for  hills  which  are  long  and  narrow  in 
the  flow  direction,  the  height  of  the  near-surface  velocity  minimum 
differs  greatly  from  that  predicted  by  the  Froude  criterion.  Since 
neither  of  the  two  velocity  minima  need  equal  zero,  a  CSH  is  not 
defined  and  layer  partitioning  cannot  proceed  in  such  a  straight 
forward  manner.  Instead,  it  appears  that  momentum  equations  which 
include  pressure  gradient  forces  may  also  be  required  to  account 
for  the  added  flow  complexities  (see  LINCOM-T  discussion) . 

Moving  on  to  other  issues,  in  CALMET  horizontal  divergence  is 
minimized  after  kinematically  adjusting  the  vertical  velocities. 
In  option  b  another  model,  usually  CSUMM,  an  old  hydrostatic 
version  of  CSU  RAMS,  provides  step  1.  In  step  2  CALMET  includes  a 
1/r^  objective  analysis  and  smooths  the  resulting  field.  Mass 
conservation  supplies  vertical  velocity,  w,  through  the  standard 
shallow  convection  form  of  the  continuity  equation,  dujdx-^  -  0  ,  or 
if  the  model  top  is  subsequently  forced  to  the  condition,  w  =  0, 
divergence  is  minimized  once  again. 

The  kinematic  and  slope  flow  adjustments  seem  quick  but  ad  hoc. 
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The  slope  flow  adjustment  requires  that  the  domain  mean  wind  be 
defined.  Studies  in  the  Los  Angeles  basin  suggest  that  when  the 
domain  mean  wind  is  difficult  to  specify,  as  is  often  the  case  for 
complex  regional  scale  terrain,  CALMET's  slope  flow  and  CSH 
adjustments  may  indicate  wind  directions  quite  opposite  those 
measured;  so  the  step  1  wind  field  at  times  does  not  even  resemble 
the  tower  winds  (GARB  and  UC  Davis,  personal  communications) . 
However,  specifying  a  domain  mean  wind  is  less  problematic  for 
smaller  domains  such  as  Vandenberg.  And  the  objective  analysis 
included  in  the  step  2  adjustment  will  re-align  the  winds  with  the 
available  tower  vectors. 

CALMET  does  use  different  approaches  for  surface  energy  balance, 
and  hence  heat  flux  and  stability,  depending  on  overwater  or 
overland  conditions.  But  at  the  same  time,  as  in  all  dry 
diagnostic  models,  the  effect  of  the  land/sea  and  cloud/clear 
interfaces  is  not  explicitly  addressed  in  its  entirety.  The  edge 
of  the  Pacific  coast  stratus  deck  hov'ers  chronically  over 
Vandenberg  much  of  the  year.  Plume  dispersion  usually  changes 
dramatically  across  the  edge  of  stratus  decks,  due  to  pressure 
induced  wind  accelerations  and  secondary  circulations,  sudden  jumps 
in  inversion  height,  and  increased  turbulence  on  the  sunny  side. 
At  the  cloud  edge  the  secondary  circulation  seems  to  involve 
vigorous  boundary  layer  scale  motions  which  may  lead  to  much  more 
rapid  vertical  mixing.  This  is  consistent  with  the  results  of  the 
Lompoc  Valley  Diffusion  Experiment  (Skupniewicz  et  al.  1992)  which 
showed  essentially  complete  vertical  mixing  within  the  boundary 
just  beyond  the  stratus  edge.  Though  stability  changes  which  drive 
the  puff  model  can  partly  account  for  some  of  these  effects,  the 
CALMET  wind  fields  are  subject  to  layer-segregated  divergence 
minimization.  So  the  model  cannot  effectively  treat  boundary  layer 
scale  vertical  miotions  encountered  across  the  stratus  edge. 
Changes  across  the  coastline  itself  are  similar  but  less  dramatic 
(Skupniewicz  et  al.  1991b). 

With  regard  to  other  issues,  the  assumed  exponential  decay  of 
vertical  velocity  in  the  kinematic  adjustment  does  not  account  for 
the  effects  of  the  strong  low  subsidence  inversions  endemic  to 
Vandenberg.  During  daytime,  tnese  inversions  tend  to  be  terrain 
parallel  more  than  flat,  resulting  in  substantial  vertical 
velocities  at  the  top  of  as  well  as  within  the  boundary  layer 
(Kamada  et  al.  1990a,  b  and  Skupniewicz  et  al.  1991b).  In  CALMET 
vertical  velocity  does  not  appear  to  be  readjusted  after  making  the 
slope  flow  and  blocking  adjustments  to  the  horizontal  winds.  This 
is  mildly  puzzling  but  perhaps  peripheral,  since  the  horizontal 
winds  are  what  are  really  desired. 

CALMET  introduces  spatio/temporal  variability  of  the  inversion 
height,  but  the  inversion-height-change-with-time  estimate  (eqn. 
2-42)  might  be  revised,  since  it  applies  to  entrainment-induced 
rather  than  subsidence  inversions.  Vandenberg 's  typical 
temperature  jumps  of  10  -  15 °C  across  the  subsidence  inversion  will 
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make  the  square  root  term  in  eqn.  2-4  2  imaginary.  Perhaps  in 
the  denominator  should  be  reinterpreted  as  the  lapse  rate  within 
rather  than  above  the  inversion  zone.  Alternatives  are  considered 
in  Kamada  (1988a, b)  and  Kamada  et  al.  (1989,  1990a, b,  1992a) .  With 
regard  to  eqn.  2-44,  Arya's  (1982)  study  of  neutral  boundary  layer 
inversion  height  formulae  conclude  that  typical  errors  are  on  the 
order  of  50  -  100%.  Another  possibility  is  Troen  and  Mahrt's 
(1986)  single  formula  which  applies  to  all  stabilities.  CALMET 
supplies  another  ad  hoc  form,  eqn.  2-45,  for  advective  changes. 
But  by  now  we  suspect  that  a  single  prognostic  inversion  height 
equation  which  includes  heating,  subsidence,  and  advection  might  be 
preferable.  An  example  implementation  is  given  in  Glendening  et 
al.  (1986). 

For  other  particulars,  in  eqn.  2-11  it  is  unclear  how  upslope  flow 
is  generated  (S  >  0),  for  negative  temperature  lapse  rates  near  the 
surface,  since  this  makes  the  square  root  imaginary.  The  surface 
energy  balance  is  based  on  Holtslag  and  van  Ulden,  (1982)  and 
supplies  CAIiMET's  initial  surface  layer  wind  and  temperature 
profile.  However,  the  user's  guide  provides  no  description  for 
obtaining  the  reqviired  incoming  and  outgoing  long-wave  radiation  in 
eqn.  2-29.  We  assume  that  this  is  also  taken  from  Holtslag  and  van 
Ulden.  More  sophisticated  radiation  schemes  are  also  presented  in 
Kamada  and  Flocchini  ( 1984a, b  and  1986)  .  Long-wave  radiation  is  a 
strong  function  of  cloud  cover  Our  own  experience  with  eqn.  2-39 
is  that  it  is  too  crude  to  be  useful  in  estimating  stability  for 
dispersion  purposes,  sometimes  leads  to  non-convergence  in  L,  the 
Obukhov  length,  and  requires  the  coefficient  to  be  tuned  for  the 
particular  site.  As  an  alternative,  we  suggest  the  energy  balance 
method  given  in  Appendix  A  of  Skupniewicz  et  al.  (1992) .  Also, 
the  Deardorff  convective  velocity  scale  in  eqn.  2-48  may  not  apply 
as  well  as  the  shear  inclusive  velocity  scale  given  in  Kamada 
(1992b)  to  Vandenberg's  usual  moderately  convective,  strongly 
baroclinic  conditions. 

RUN-TIME: 

We  suspect  that  most  of  CALMET 's  rather  lengthy  time  is  spent  with 
objective  analysis  and  iterative  divergence  minimization.  We 
cannot  determine  from  the  description  whether  the  divergence 
minimization  is  fully  three-dimensional  (see  MATHEW  review  for 
equations) ,  or  uses  the  faster  pseudo-two-dimensional  equation  over 
several  flow  surfaces,  such  as  in  WOCSS  and  MACHWIND  (reviewed 
below) .  Apparently  run-time  is  also  a  strong  function  of  the 
number  of  surface  stations  involved  in  the  analysis  (CARB,  personal 
communication,  1992).  For  a  30x30x5  (4,500  grid  point)  domain  at 
4  kilometer  resolution  using  14  surface  stations,  UC  Davis  reports 
about  ten  minutes  for  execution  time  on  a  VAX  3200  workstation,  a 
"12  MIPS  machine.  Suppose  three  of  the  five  minutes  allotted  for 
THC  output  is  available  to  compute  the  windfield.  Vandenberg's 
requirements  are  for  a  50km  x  50km  domain  at  500m  resolution  with 


20 


at  least  four  levels  (40,000  grid  points)  incorporating  30  towers 
and  five  SODAR.  For  this  configuration  we  expect  emergency 
response  using  CALMET  to  require  a  400  MIPS  computer  (perhaps 
available  in  3  years  for  under  $50K) .  Currently  available 
workstations  priced  less  than  $50K  run  at  speeds  of  50  -  120  MIPS 
and  perhaps  5-30  MFLOPS  (million  floating  point  instructions  per 
second)  as  fig.  1  indicates.  Thus,  for  now  it  appears  that 
CALMET /CALPUFF  is  really  a  regional  air  quality  model  intended  for 
“200km  X  300km  domains  with  a  4  km  mesh,  rather  than  a  high 
resolution  local  scale  emergency  response  model.  Some  sub-grid 
scale  terrain  effects  can  be  included  in  CALPUFF,  rather  than 
CALMET,  by  modeling  sub-grid  scale  hills  as  ellipsoids.  This  may 
be  useful  for  regional  air  quality  studies,  but  we  question 
whether,  in  lieu  of  500  meter  resolution,  this  approach  is  adequate 
for  Vandenberg. 

For  one  thing,  4  km  resolution  of  the  coastal  interface  (either 
land/sea  or  cloudy/clear)  is  probably  inadequate.  Much  of  the 
Vandenberg  coastline  terrain  consists  of  abrupt  bluffs  rather  than 
smooth  beach;  so  it  is  kinematically  and  mechanically  quite 
disruptive  of  the  near  surface  flow.  The  Hypergolic  Stockpile  and 
Storage  Facility  (HSSF)  and  other  potential  "cold  spill"  sites  lie 
within  two  kilometers  of  the  coastline.  For  roughly  6  months  of 
the  year  the  seabreeze  blows  the  stratus  deck  inland,  while  solar 
heating  burns  it  back  to  within  a  few  kilometers  of  the  coast.  So 
for  local  scale  dispersion,  these  large  near-coastal  changes 
probably  compel  1  km  accuracy  in  locating  the  coastal  and  stratus 
edge  interfaces,  as  well  as  paramtrizations  specific  to  the  cloud 
edge.  Skupniewicz  et  al.  (1992)  describe  such  a  two-zone 
convective  scaling  procedure  and  algorithms  for  high  resolution 
specification  of  the  cloudy/clear  interface  from  GOES  and  AVHRR 
satellite  images. 

CALPUFF: 

INTRO: 

CALPUFF  seems  to  be  a  sophisticated,  fairly  inclusive,  puff  model 
which  features  multiple  puffs  and  alterations  due  to  building 
down-wash  near  the  source,  transitional  plume  rise,  vertical  wind 
shear  and  subgrid-scale  terrain  interactions.  It  also  includes  wet 
and  dry  deposition,  .'ir  quality  oriented  chemistry,  overwater 
transport,  coastal  interaction,  some  options  for  dispersion 
coefficients,  and  near-field  puff  stretching.  Puff  splitting  is 
also  included  to  allow  for  wind  shear  and  plume  bifurcation  in 
complex  terrain. 

PHYSICS: 

CALPUFF's  near-field  puff  elongation  (slug  mode)  probably  accrues 
large  savings  in  CPU  time  for  commensurate  accuracy.  That  is,  for 
standard  puff  models  small  puffs  (<  100  meter  radius)  must  be 
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released  frequently  to  maintain  near-field  accuracy  in  a  time- 
varying  windfield.  So  to  avoid  artificial  concentration 
oscillations,  the  mean  distance  between  puffs  at  the  time  of 
release  should  not  exceed  1.5  puff  radii. 

CALPUFF  has  three  options  for  dispersion  coefficients:  direct 
measurement  of  turbulence  intensities,  similarity  expressions,  and 
Pasquill-Gif ford  stability  classes.  Obtaining  accurate  direct 
measurements  of  vertical  turbulence  intensity  is  unlikely,  so 
CALPUFF 's  default  option  is  to  use  similarity  based  expressions. 
These  are  taken  from  Panofsky  (1977),  Hicks  (1985),  Arya  (1984), 
Nieuwstadt  (1984),  Deardorff  (1975),  and  the  empirical 
interpolation  formulae  of  Scire  et  al.  (1990)  .  Analyzing  these 
options,  we  rote  that  data  under  unstable  conditions  show  that 
plumes  from  surface  releases  tend  to  lift  rather  quickly  above 
ground  in  the  near-field,  leaving  ground  level  concentrations  much 
lower  than  is  indicated  by  parametrizations  using  Pasquill-Gif ford 
stability  classes  A  and  B  (Willis  and  Deardorff,  1976,  1978;  Lamb, 
1978,  1979;  Nieuwstadt,  1980;  and  Briggs,  1985,  1986).  On  the 
other  hand,  the  congruence  between  data  and  Nieuwstadt 's  similarity 
expressions  for  stable  cases  has  been  questioned;  the  alternative 
scaling  by  Sorbjan  (1989)  may  be  more  robust.  For  unstable  to 
neutral  cases,  an  alternative  to  empirical  interpolation  is  the 
Kamada  (1992b)  scaling  which  1)  is  physically  derived  and  2)  also 
treats  mechanical  turbulence  due  to  entrainment,  baroclinicity ,  and 
surface  layer  shear. 

Another  issue  is  that  CALPUFF  lacks  spectral  filtering.  Without 
spectral  filtering,  a  puff  model  designed  for  use  with  a  static 
wind  field  will  over-estimate  dispersion,  if  that  wind  field  is 
frequently  updated.  The  point  is  that  eddies  much  larger  than  the 
puff  dimension  will  induce  meander  rather  than  puff  growth.  If 
time-averaged  which  include  meander  are  used,  then  diffusion 
will  be  exaggerated.  Since  THCs  must  account  for  plume  meander 
during  the  dispersion  period,  spectral  filtering  is  not  needed  with 
fixed  winds,  provided  the  Cy ^  averaging  times  and  hazard  dispersion 
periods  are  about  equal.  But  for  CSUMM's  time-varying  winds  or  for 
frequently  updated  met  data,  meander  must  be  filtered  from  the  puff 
model  to  avoid  double  counting  it  in  both  the  windfield  and  puff. 
The  band  width  for  the  required  high-pass  spectral  filtering  must 
also  change  with  individual  puff  size.  Stability  and  wind 
direction  dependent  velocity  spectra  for  Vandenberg  were  presented 
for  this  purpose  in  Skupniewicz  et  al.  (1989) . 

An  issue  already  addressed  in  the  CALMET  review  is  that  within  one 
grid  cell,  the  transition  from  marine  to  continental  dispersion 
rates  can  be  complicated  by  the  presence  of  both  the  coastline,  as 
well  as  cloud/clear  interfaces  (Skupniewicz  et  al.  1991b).  Since 
much  of  south  Vandenberg 's  hilly  terrain  simply  extends  into  the 
sea,  the  coastline  discontinuuity  is  quite  disruptive  of  the  near 
surface  flow.  Not  only  is  there  mean  flow  acceleration  but  also 
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augmented  turbulence  due  to  large  velocity  gradients  and 
significant  turbulence  anisotropy,  all  initiated  "  1  -  2  km  upwind 
of  the  potential  release  sites.  This  may  occur  together  with  even 
more  dramatic  flow  disruptions  across  a  stratus  edge  which  often 
lies  at  or  just  beyond  the  release  site.  The  stability  indices  and 
sub-grid  procedures  in  CALPUFF  must  account  to  some  degree  for 
these  features,  but  again  this  seems  to  demand  high  resolution  and 
specific  validation  for  Vandenberg.  Current  data  is  most  valuable, 
so  maintenance  of  tower  057  next  to  the  HSSF  and  the  SODAR  and 
tower  at  Building  900  is  critical  for  HSSF  releases.  But  bear  in 
mind  that  the  horizontal  dispersion  coefficients,  and  a^, 
obtained  from  SODARS  have  been  considered  unreliable  (Neff,  1986) . 

Again,  CALPUFF 's  CSH  algorithm  is  consistent  with  theories  of 
stratified  flow  by  Hunt  (1978)  and  Carruthers  and  Hunt  (1990)  . 
These  have  achieved  some  empirical  but  not  theoretical  confirmation 
(Smith,  1990) .  Allowance  for  entrainment  and  internal  boundary 
layer  growth  implies  that  downwind  fumigation  can  be  treated.  The 
chemistry  model  is  interesting  but  not  quite  relevant,  since  the 
pollutants  of  general  air  quality  concern  are  different  from  the 
propellants  and  oxidizers  used  near  the  launch  complexes.  Relevant 
chemistry  must  be  installed  in  CALPUFF,  perhaps  from  the  REED-M  and 
AFTOX  chemistry  modules.  On  the  other  hand,  many  of  the  reviewed 
diffusion  models  lack  any  treatment  of  downwind  chemistry.  So 
CALPUFF  may  have  some  advantage  in  this  regard. 

CALMET/ CALPUFF  SUMMARY: 

Apart  from  CALPUlF's  technical  merit  in  the  context  of  a  static 
wind  field,  a  major  strength  of  CALMET /CALPUFF  seems  to  be  that  it 
is  a  fairly  comprehensive,  modular,  and  operational  system 
maintained  by  an  interested  public  agency,  the  California  Air 
Resources  Board  (CARS) .  The  entire  code  and  documentation  is 
available  free  from  EPA  electronic  bulletin  boards  or  through  CARB. 
and  is  undergoing  an  extensive  review  for  regional  air  quality 
purposes  at  the  University  of  California  at  Davis.  Apparently,  no 
site  license  is  required,  nor  are  there  any  special  restrictions  on 
local  modification.  On-line  technical  support  may  be  available. 
However,  this  issue  alv/ays  awaits  de  facto  confirmation.  Since  the 
code  is  actively  supported  by  a  large  local  agency  and  widely 
available,  the  climate  for  general  acceptance  may  be  good,  time  to 
obsolescence  relatively  long,  and  life-cycle  costs  commensurate ly 
low.  CARB  also  plans  validation  studies  using  its  Los  Angeles  Air 
Basin  (1987)  and  San  Joaquin  Valley  (1990)  tracer  data  sets  (CARB, 
personal  communication) .  Such  studies  should  be  regarded  as 
supplementary  to,  rather  than  a  substitute  for  studies  specific  to 
Vandenberg,  since  the  domain  scale,  station  densities,  and  terrain 
differ  enough  to  effect  qualitative  as  well  as  quantitative  changes 
in  model  performance.  Suitably  tailored  for  Vandenberg 's  needs, 
the  cost,  time,  and  training  involved  in  procurement,  testing, 
validation,  and  operational  installation  should  be  similar  to  that 
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of  other  models  not  specifically  tailored  for  the  MARSS  system. 
For  most  models,  installation  of  Vandenberg  data  bases  and 
appropriate  sizing  and  resolution  of  the  grid  is  likely  to  be 
considerably  more  v.ime  consuming  than  the  relatively  minor  effort 
involved  in  porting  a  Fortran  based  model  to  MARSS. 

Compared  to  an  ideal  Vandenberg  model,  chief  weaknesses  and 
omissions  seem  to  be  near-field  dense  gas  effects,  liquid  spill 
source  term  assessment,  CALMET's  slow  speed,  inadequate  grid 
resolution,  ad  hoc  physics,  incomplete  coastal  and  cloud/clear 
interface  parametrizations,  and  lack  of  Vandenberg  tracer 
validation.  In  particular  CALMET's  slowness  makes  it  unsuitable 
for  Vandenberg  emergency  response.  However,  for  the  step  1  CALMET 
wind  field,  option  b  allows  substitution  with  another  wind  field 
model.  If  CALPUFF  is  chosen,  we  recommend  substituting  CALMET  with 
a  faster,  higher  resolution  model,  such  as  LINCOM  and  MACHWIND.  If 
CALPUFF  is  installed,  we  also  suggest  that  ir  be  retrofitted  to  use 
filtered  velocity  spectra  for  its  dispersion  coefficients.  This 
may  be  a  difficult  task  requiring  fundamental  changes  to  puff 
model's  physical  assumptions  as  well  as  coding.  Perhaps  the  AFTOX 
or  REED-M  source  term  modules  could  also  be  modified  for  inclusion 
in  CALPUFF  or  any  of  the  other  reviewed  models,  since  they  all  seem 
to  lack  this  feature. 
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PGEMS,  MELSAR/MESOI  (Allwine,  Rainsdell,  et  al.): 

INTRO: 

The  PGEMS  model  is  slated  to  become  operational  at  the  Diablo 
Canyon  nuclear  power  plant.  The  core  of  the  modeling  system 
consists  of  the  wind  field  portion  of  MELSAR  and  the  MESOI  puff 
diffusion  model  (both  developed  by  Battelle  Labs) .  The  Diablo 
Canyon  domain  is  quite  similar  to  south  Vandenberg ' s .  Prominent 
features  are  a  turning  coastline,  rugged  coastal  hills  with  peaks 
up  to  ~  400  meters,  a  number  of  intervening  canyons,  and  a  broad 
series  of  valleys  surrounding  the  coastal  hills  from  the  northwest 
Morro  Bay  region  to  the  southeast  Avila  Beach  region.  Since  a 
validation  study  was  undertaken  for  this  area,  one  initial  hope  was 
to  avoid  the  expense  of  a  separate  study  for  Vandenberg,  should 
PGEMS  be  seriously  considered.  We  discuss  this  possibility  below. 

MELSAR  wind  field  (K.J.  Allwine  and  C.D.  Whiteman): 

PHYSICS: 

Following  a  Cartesian  to  terrain  following  grid  transformation, 
mass  consistency  is  based  on  the  standard  shallow  convection  form 
of  the  continuity  equation.  Some  explanation  seems  missing  for 
eqns.  2-12  to  2-16  of  the  MELSAR  technical  documentation  (PNL-5460 
Vol.l  UC-1) ,  since  units  are  not  in  agreement  for  the  left  and 
right  hand  sides.  These  equations  dictate  grid  transforms  for  the 
wind  field.  Nine  vertical  levels  are  specifiable.  However,  we 
cannot  reasonably  interpret  the  sentence  on  p.  2-20,  "These  input 
upper-air  observations  at  each  gamma  level  are  then  spatially 
interpolated  using  a  1/r*  weighting  to  each  surface  weather 
station,  with  the  ground  level  interpolated  winds  being  replaced  by 
the  surface-wind  observations".  The  remainder  of  this  part  of  the 
description  is  equally  murky,  but  a  reasonable  interpretation  would 
be  that  the  objectively  analyzed  vertical  v'inds  are  held  fixed, 
while  the  continuity  equation  is  used  to  adjust  the  horizontal 
winds  to  he  divergence  free.  The  wind  field  for  each  surface  is 
represented  as  sums  of  amplitude  functions  multiplied  by  fifth 
degree  orthonormal  Chebyshev  polynomial  basis  functions,  according 
to  an  undescribed  non-linear  least  squares  technique.  The  purpose 
is  to  greatly  reduce  memory  storage  requirements.  The  solution 
technique  assumes  that  there  are  no  domain  scale  horizontal  wind 
gradients  and  no  general  subsidence  or  lifting.  These  are  not  good 
assumptions  for  Vandenberg.  Again,  layer-segregated  mass 
continuity  does  not  adequately  describe  non-terrain  induced, 
boundary  layer  scale  vertical  motions  in  unstable  atmospheres, 
convective  clouds,  or  at  the  stratus  cloud  edges  which  are  endemic 
to  the  Vandenberg  domain. 

A  CSH  is  computed  to  determine  terrain  induced  flow  distortion 
under  stable  conditions,  but  details  are  puzzling.  Temperature 
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lapse  rates  for  the  hill  Froude  number,  s  UqH/N  <  1  (see  CALMET 
review  for  more  detail) ,  are  determined  from  input  temperature/ 
pressure  soundings  and  surface  observations.  As  shown  by  Hanna 
(1987)  and  specifically  for  Vandenherg  by  Skupniewic?;  et  al. 
(1992) ,  energy  balance  methods  are  more  accurate  than  actual 
surface  observations  in  estimating  surface  layer  lapse  rates. 
Though  the  treatment  of  temporal  variation  in  lapse  rate  seems 
reasonable,  the  description  of  averaged  surface  and  obstacle  height 
lapse  rates  from  soundings  is  not  clear  or  convincing. 

The  mixing  height  is  assumed  to  be  the  maximum  of  the  mechanical 
and  convective  mixing  heights.  The  mechanical  mixing  height  is 
given  as  0.0053  Ug,  where  Ug  is  the  free  stream  velocity  (m/s). 
Obviously,  the  constant  is  in  error.  Convective  mixing  heights  are 
estimated  by  drawing  a  vertical  line  from  the  current  surface 
potential  temperature  up  to  where  it  intersects  the  morning 
sounding.  So  this  assumes  that  surface  heating  variations  perturb 
the  otherwise  flat  boundary  layer  top.  Lapse  rates  and  mixing 
heights  are  interpolated  from  sounding  and  surface  station  sites  to 
a  10  X  10  grid  spanning  the  domain.  Upon  reaching  the  strong 
subsidence  inversion,  this  method  implies  that  boundary  layer 
heights  become  more  uniform,  i.e.,  less  than  terrain  following, 
contrary  to  the  day-time  studies  of  Kamada,  et  al.  (1990a, b). 

Pasquill-Gif ford  stability  class  is  estimated  for  each  gridpoint 
and  hour  according  to  wind  speed  and  convective  mixing  height.  A 
thermodynamic  method  is  presented  which  includes  a  topographic 
amplification  factor  (TAF)  to  estimate  coupling  between  valley  and 
ambient  flows.  Equations,  2-23  to  2-26  prognose  the  convective 
boundary  layer  height  and  inversion  top  height.  Together  these 
define  the  inversion  depth  within  a  valley.  They  are  in  turn  used 
to  define  an  hourly  flow  coupling  coefficient.  At  Vandenberg's 
required  fine  grid  resolutions,  this  may  not  be  needed  for  the 
Lompoc  Valley.  However,  for  the  given  inland  heating  and  distance 
from  the  coastline  or  stratus  cloud  edge,  Lompoc  measur'^ments  show 
inversions  much  lower  than  expected  (Skupniewicz  et  al.,  1990, 
1992)  .  This  may  stem  from  1)  an  inversion  deepened  by  cold 
nocturnal  drainage  into  Lompoc  Valley,  and  2)  subsidence  in  the 
valley  center  which  compensates  for  day-time  upslope  flow  along  the 
sides,  effects  included  in  eqns.  2-23  to  2-26.  The  method  is  based 
on  studies  involving  topographical  features  much  larger  than  south 
Vandenberg's  canyons.  Since  pressure  gradient  driven  slope  flows 
increase  as  the  square  of  the  horizontal  scale  over  which  the 
gradient  applies,  flows  over  short  slopes  may  be  limited  (Mikkelsen 
and  Kamada,  unpublished).  This  would  reduce  the  method's 
quantitative  accuracy  for  Honda  and  other  local  canyons.  Yet  a 
number  of  cases  during  the  Mt.  Iron  tracer  study  showed  flow 
decoupling  over  Honda  Canyon  (Kamada  et  al.  1992c, d) .  So  the 
MELSAR  method  invites  further  scrutiny  for  both  the  Lompoc  Valley 
case  and  sub-grid  scale  canyon  flows. 
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The  physical  basis  of.  the  MELSAR  wind  model  is  much  simpler  than 
CALMET,  LINCOM,  or  even  WOCSS.  In  the  context  of  computer  power 
available  or  soon  to  be  available  at  Vandenberg,  the  question  is 
whether  the  presumed  speed  advantage  offsets  the  presumed  accuracy 
loss.  Tower/ SODAR  data  density  and  the  quality  of  objective 
analysis  are  more  critical  in  assuring  the  accuracy  of  output  wind 
fields  when  using  less  physically  based  models  such  as  MELSAR. 
Improved  objective  analyses  are  available  as  in  LINCOM.  But  bear 
in  mind  that  most  of  Vandenberg *s  many  met  sensors  are  sited  near 
the  coast.  The  NOAA  buoys,  oil  platforms,  and  Santa  Barbara  Air 
Quality  Management  towers  are  not  included  in  the  on-line  data 
stream.  We  have  not  found  discussions  of  run-time  in  the  available 
PGEMS  documentation  for  the  complete  model.  Glantz  and  Burk  (1990) 
mention  that  computation  of  just  the  wind  field  alone  requires  "  30 
seconds  per  hour  of  simulation  on  a  MicroVax  II.  Porting  MELSAR  to 
a  current  486  PC  might  reduce  run-time  to  less  than  one  second. 
However,  for  Vandenberg 's  purposes,  enhanced  accuracy  is  probably 
more  important  than  run-times  much  below  30  seconds. 

VALIDATION: 

A  wind  model  verification  project  is  described  in  Appendix  E  of 
Phase  I  of  the  PGEMS  technical  document,  referred  to  as  Appendix  A 
of  Volume  I  of  Phase  III  of  the  PG&E  Report  009.5-88.4  (Thullier  et 
al.  1988).  Seventeen  hourly  wind  fields  were  used  from  nine 
surface  stations  and  two  SODAR  sites  for  5:00  to  21:00  June  11, 
1985.  The  approach  was  to  compare  wind  fields  with/without  data 
from  a  given  site  and  with/without  the  critical  streamline  height 
(CSH)  option.  Only  two  stations  were  tested,  the  San  Luis  Obispo 
APC  site  and  VC  Doppler  station.  Results  for  the  six  runs  were 
decidedly  mixed  with  disagreement  between  modeled  and  measured 
winds  ranging  from  5°  and  0.3m/s  to  41°  and  l.lm/s.  Use  of  the  CSH 
alternately  improved  and  degraded  agreement.  Separation  distance 
between  stations  was  rather  large  by  Vandenberg  standards,  on  the 
order  of  a  ten  to  twenty  kilometerss.  For  Vandenberg  this  test  was 
too  cursory  to  validate  the  wind  flow  model.  Little  was  said  about 
the  local  terrain  surrounding  the  two  station  sites  in  question,  so 
it  was  not  possible  to  evaluate  the  strength  of  the  flow  model  on 
this  basis.  Congruent  with  these  thoughts  the  evaluation  portion 
of  the  report  (Phase  III,  pps .  24  -  27)  suggests  that  the  wind 
field  relies  strongly  on  interpolations  from  individual  site 
measurements  and  notes  that  "mass  consistent  adjustments  are  of 
limited  value  in  the  absence  of  measurements"  and  that  "the 
diagnostic  model  is  simply  an  interpolative  scheme  that  attempts  to 
flesh  out  the  details  of  the  wind  field  in  a  manner  consistent  with 
the  principle  of  mass  continuity"  (p.  62).  A  look  at  fig.  8,  p.  25 
of  the  report  suggests  that  without  the  CSH  option  the  model  will 
ignore  most  terrain  effects  beyond  that  objectively  analyzed  from 
station  data.  So  high  tower  and  SODAR  density  are  critical  to 
MELSAR 's  use  in  complex  terrain.  The  absence  of  Lowers  along  the 
northwestern  promontory  and  in  See  Canyon  meant  that  terrain 
channeling  and  strong  upslope  flow  in  those  areas  were  not 
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predicted.  Thus,  major  discrepancies  appeared  in  the  dispersion 
predictions  on  some  occasions.  On  the  other  hand,  towers  may  be 
aliased  by  local  terrain,  resulting  in  local  flows  not  indicative 
of  larger  scales.  This  was  seen  in  the  Mt.  Iron  vs.  LINCOM/ 
RIMPUFF  results  for  flow  over  Vandenberg's  Honda  Canyon  (Kamada  et 
al.  1992c, d) .  This  highlights  the  need  for  inclusion  of  as  much 
physics  in  wind  flow  models  as  possible,  within  the  available  time 
and  computer  power  constraints. 

A  more  extensive  tracer  study  took  place  from  the  power  plant  and 
Morro  Bay  for  seven  non-consecutive  days  during  the  period  August 
31  to  September  15,  1986.  Releases  occurred  from  0800  to  1600  LDT 
and  sampling  from  0700  to  1900  LDT.  This  sampling  period  and 
schedule  is  quite  comparable  to  the  1989  Lompoc  Valley  Diffusion 
Study  for  releases  from  the  HSSF  and  the  mouth  of  the  Lompoc  Valley 
at  south  Vandenberg.  However,  the  range  of  wind  patterns  from  the 
power  plant  were  considerably  more  diverse  than  seen  from  the  HSSF. 
The  two  modeled  Diablo  Canyon  domains,  northwesterly  and 
southeasterly,  were  both  50  x  50  km  at  1  km  grid  spacing  and  were 
offset  from  each  other  by  10  km  along  the  north/south  axis  and  12 
km  in  the  east/west.  The  northwesterly  domain  was  used  to  model 
southeasterly  flow  from  the  releases  at  the  power  plant,  while  the 
southeasterly  domain  was  used  for  northwesterly  flows.  Three 
SODAR,  eighteen  surface  stations,  and  about  150  bag  samplers  were 
employed  in  the  study.  The  CSH  option  was  not  employed  for  the 
windfield  model  runs.  Twelve  puffs  were  released  per  hour.  The 
model  computes  concentrations  on  a  31  x  30  grid  encompassing  the 
domain  and  at  up  to  25  receptor  sites,  so  eleven  separate  modeling 
runs  were  employed  for  each  release  to  treat  the  150  receptors.  A 
terrain  parallel  upper  boundary  at  1200  m  AGL  was  employed.  Wind 
fields  were  updated  at  15  minute  intervals.  In  an  updated  PGEMS  v. 
1.1  the  list  has  been  expanded  to  include  250  receptor  sites,  so 
separate  runs  are  no  longer  necessary. 

Graphical  results  do  not  immediately  appear  as  congruent  as  the 
Vandenberg  Mt.  Iron  vs.  LINCOM/RIMPUFF  set.  But  this  may  be 
misleading  because  Diablo  Canyon  involved  much  longer  ranges,  along 
with  a  lower  relative  density  of  wind  measurements.  In  contrast  to 
the  rather  cursory  wind  field  analysis,  the  report  gives  a  lengthy 
statistical  and  graphical  analysis  to  show  the  considerable 
disparity  between  near-field  modeled  and  measured  concentrations. 
The  report  assumes  that  building  down-wash  and  misleading  wind 
direction  fluctuation  statistics  (sigma  0)  from  the  release  site 
tower  lead  to  more  stable  PG  diffusion  classes  than  were  effective 
downwind  or  above  the  surface  layer.  Apparently,  a  building  down- 
wash  algorithm  has  since  been  installed  in  the  puff  model.  However, 
as  noted  above,  the  remaining  issue  is  that  surface  releases  tend 
to  loft  quickly  in  unstable  conditions.  So  it  is  well-established 
that  PG  classes  A  and  B  will  overpredict  near-field  ground  level 
doses.  We  recommend  deriving  dispersion  coefficients  from  surface 
energy  balance  and  similarity  theory,  rather  than  tower  wind/ 
temperature  profiles  such  as  in  PGEMS,  if  l  ability  and  wind 
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dependent  velocity  spectra  are  unavailable  at  the  release  site  and 
for  each  tower. 

We  are  puzzled  that  results  from  some  runs  using  the  CSH  option 
were  not  presented.  However,  roost  of  the  releases  occurred  for 
atmospheric  stabilities  too  weak  to  warrant  its  use.  The  report 
makes  the  general  comment  that  the  windfield  agreement  with  data 
tended  to  degrade  near  the  domain  boundaries. 

Far-field  results  were  rather  poor,  but  comparable  to  other  complex 
terrain  studies.  "Correlations  between  predicted  and  observed 
values  were  generally  insignificant  (p.74)".  A  distinct  bias 
toward  overprediction  of  concentrations  appeared,  perhaps  partly 
due  to  the  lack  of  filtered  wind  velocity  spectra  for  dispersion 
coefficients  (corrected  in  v.  l.l).  MESOI's  lack  of  puff  splitting 
or  much  physics  in  the  wind  field  may  be  partly  responsible  also. 
Vertical  and  horizontal  puff  splitting  were  also  added  to  PGEMS  v. 
1.1.  So  improvements  have  been  made,  but  there  remains  some 
necessary  disagreement  between  predicted  and  observed  results 
stemming  from  the  stochastic  nature  of  atmospheric  turbulence.  So 
a  high  degree  of  far-field  correlation  cannot  be  expected, 
particularly  in  complex  terrain.  However,  the  Diablo  Canyon  study 
results  do  not  warrant  acceptance  of  the  MELSAR  wind  flow  model  for 
operational  use  at  Vandenberg  without  additional  on-site  studies. 


MESOI  PUFF  MODEL  (J.V.  Ramsdell,  G.F.  Athey,  C.S.  Glantz) ; 
INTRO : 

Information  for  this  review  comes  from  an  NRC  report  published  by 
the  above  authors  in  1983.  We  suspect  that  further  work  must  have 
proceeded  on  this  model  in  the  interim.  Installation  of  certain 
improvements  were  surmised  from  comments  in  other  MELSAR 
documentation.  There  perhaps  are  others  which  we  have  not 
detected.  Be  that  as  it  may,  from  our  best  information  MESOI  seems 
to  be  a  standard  puff  diffusion  model  which  fairly  mature  in  the 
sense  that  itt  also  incorporates  puff-splitting,  plume  rise, 
building  down-wash,  wet  and  dry  deposition,  and  half-life  decay 
rates  for  radioactive  materials  and  chemicals. 


PHYSICS: 

Plume  rise  is  computed  from  the  standard  effective  stack  height 
equations  of  Briggs.  These  equations  are  not  recommended  if  the 
release  does  not  occur  through  a  stack.  For  puff  centers  the 
resulting  effective  release  height  parallels  the  terrain.  Away 
from  puff  centers  the  effective  release  height  remains  fixed  or 
matches  terrain  height,  if  the  terrain  height  happens  to  exceed  it. 
Terrain  height  is  computed  by  interpolating  gridded  terrain  data. 
The  inversion  is  set  to  be  terrain  parallel.  Total  reflection  from 
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the  inversion  is  assumed,  so  partial  puff  penetration  is  not 
allowed.  However,  if  the  release  occurs  above  the  mixing  layer, 
downward  penetration  is  allowed. 

Diffusion  coefficients  are  based  on  several  data  sets  obtained  in 
the  1960 's  and  categorized  by  Pasquill  stability  classes.  We 
repeat  that  conservative  tracers  were  not  available  before  the 
advent  of  SFg  "  1970  and  that  the  A  and  B  Pasquill  stability 
classes  do  not  account  for  near-field  plume  lofting  associated  with 
convective  conditions.  So  these  features  need  updating,  if  this 
has  not  already  been  done. 

A  second  option  for  diffusion  coefficients  in  MESOI  is  to  use 
current  local  turbulence  data,  is  derived  from,  ct^/U,  the  ratio 
of  the  standard  deviation  of  cross  wind  turbulence,  to  mean  wind 
speed,  is  obtained  from  the  standard  deviation  of  the  wind 
elevation  angle  (in  radians) .  This  is  probably  a  better  method  for 
horizontal  turbulence.  However,  for  vertical  turbulence  we  suspect 
that  (j^,  if  derived  from  bivanes  not  sonic  anemometers,  probably 
does  not  reliably  measure  vertical  wind  fluctuations.  Of  course 
Vandenberg  towers  cannot  normally  be  outfitted  with  expensive 
sonics.  As  in  the  CALMET/CALPUFF  reviews,  we  again  suggest  as  a 
practical  matter  that  vertical  turbulence  be  estimated  from  surface 
energy  balance  computations  and  boundary  layer  similarity  theory, 
perhaps  in  conjunction  with  solar  and  net  radiometer  readings. 

A  technique  of  computing  virtual  source  distances  is  used  to  avoid 
sudden  changes  in  puff  dimensions  with  stability  class.  The 
spectral  filtering  caveat,  described  in  the  above  CALPUFF  review 
also  applies,  so  MESOI  may  over-estimate  diffusion,  if  the  wind 
field  is  frequently  updated  and  puff  meander  is  double-counted. 

Downwind  chemistry  ma>  be  included  within  MESOI  with  assumed  half- 
lives  for  species.  MESOI  can  also  account  for  the  creation  of 
toxic  intermediate  species,  as  well  as  final  inert  species  in  a  two 
stage  process  which  assumes  exponential  decay  and  creation  rates, 
ala  radioactive  processes.  For  the  puff  level  of  modeling,  this 
seems  appropriate.  Complete  computation  of  downwind  chemistry 
would  be  computationally  prohibitive  in  an  emergency  response 
context. 

Dry  deposition  is  estimated  using  the  standard  source  depletion 
method,  assuming  that  the  deposition  velocity  is  0.01  ms* 
regardless  of  position  or  material.  Wet  deposition  is  treated  as 
washout,  wherein  mass  loss  from  the  puffs  is  proportional  to 
precipitation  rate  and  in-puff  mass  concentrations.  Such  estimates 
are  probably  useful  to  within  an  order  of  magnitude. 

PGEMS  SUMMARY: 

In  summary  we  are  left  with  the  impression  that  MELSAR  was 
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developed  some  years  earlier  than  CALMET/CALPUFF  for  efficient  use 
on  computers  from  that  era.  Inexpensive  computer  power  has 
increased  by  two  to  three  orders  of  magnitude  in  the  interim. 
Hence,  some  of  the  simplifications  used,  particularly  in  the  MELSAR 
windfield  are  no  longer  necessary,  nor  are  they  all  appropriate  for 
Vandenberg.  MESOI  appears  to  be  a  viable  puff  model  which  includes 
computationally  simple  downwind  chemistry  capability  based  on  a 
half-life  approach,  as  well  as  many  of  the  add-ons  necessary  for 
practical  dispersion  estimates.  Unfortunately,  the  chemical  base 
does  not  include  the  types  of  rocket  propellents  of  interest  to 
Vandenberg.  As  with  CALPUFF,  the  lack  of  spectral  filtering 
registers  a  caveat,  if  MESOI  is  piggy-backed  to  a  frequently 
updated  wind  field,  as  is  expected  at  Vandenberg. 

Also,  we  find  that  errors  in  the  MELSAR  windfield  documentation 
available  to  us  have  not  been  updated  consistently-  We  do  not 
understand  some  portions  of  the  wind  field  model  description.  The 
valley  flow  theory  developed  by  Battelle  Laboratory  personnel  is 
excellent.  However,  it  is  meant  to  apply  for  grid  resolutions  that 
are  too  coarse  for  Vandenberg  use.  Again,  other  weaknesses  .are  the 
lack  of  accounting  for  near-field  dense  gas  effects  and  source 
terms  for  liquid  spills. 

On  the  other  hand  we  expect  that  MELSAR  is  already  adapted  to  DEC 
MicroVax  II  computers  using  the  VMS  operating  system,  and  Techtronix 
graphics  display  terminals;  it  is  installed  in  this  configuration 
at  the  Diablo  Canyon  nuclear  power  plant.  We  understand  that  it 
and  WCCSS  also  have  been  validated  for  Diablo  Canyon  in  similar 
coastal  terrain  less  than  100  km  north  of  Vetndenberg  and  that  the 
delivered  Battelle  model  featured  a  domain  which  could  be  varied 
between  ten  and  several  hundred  kilometers  on  a  side.  If  so,  we 
suspect  that  the  MELSAR  documentation  available  to  us  may  be 
somewhat  out  of  date,  while  an  improved  version  of  the  model  may 
exist.  Perhaps  further  contact  with  Battelle  Labs  could  clarify 
this  issue.  At  the  same  time  the  downwind  range  of  nuclear  plume 
THCs  certainly  exceeds  that  of  any  conceivable  chemical  plume 
emitted  from  Vandenberg.  So  the  emphasis  on  far-field  dispersion 
and  coarse  grids  may  be  entirely  appropriate  for  Diablo  Canyon  but 
not  for  Vandenberg.  If  so,  we  suggest  that  while  the  code  may  be 
as  accessible  as  CALMET/CALPUFF  and  the  installation  effort 
minimal,  so  might  obsolescence  time  for  such  a  model.  So  the 
overall  life  cycle  cost  of  this  modeling  system  may  be  high,  even 
if  the  initial  installation  costs  are  deceptively  low. 
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NUATMOS/CITPUFF  (D.G.  Ross  and  D.G.  Fox) 

NUATMOS  (D.G.  Ross  et  al.): 

NUATMOS  employs  terrain  following  coordinates  and  variable  vertical 
grid  spacing.  Like  CALMET,  MELSAR,  WOCSS,  MACHWIND,  and  MATHEW,  it 
is  a  mass  conservative  windflow  model.  It  differs  from  these  other 
in  that  it  seeks  to  include  dynamic  speed-up  and  stability  effects 
by  apportioning  the  adjustment  between  horizontal  and  vertical  flow 
during  iterative  divergence  minimization.  In  this  process  the 
error  functional, 
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is  minimized  by  the  variational  method,  subject  to  the  standard 
shallow  convection  continuity  constraint  (see  MATHEW  review  for 
general  description) .  Energy  conservation  arguments  and  laboratory 
data  suggested  to  the  authors  of  NUATMOS  that  the  tuning  parameter 
for  vertical  velocity  adjustments,  a'^,  is  well  characterized  by  the 
form. 
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Here  S  is  the  amount  of  speedup  over  terrain  beyond  the  initial 
upstream  velocity,  Uq,  and  Fr  =  Uq/NH  is  the  hill  Froude  number 
based  on  the  Brunt  Vaisala  frequency,  N  =  [  (g/9)  (d&/ dz)  ,  and  the 
hill  height,  H  (Ross  et  al.  1991)  (See  CALMET  and  MATHEW  reviews 
for  more  detail) .  Earlier  test  forms  for  a,  based  on  physical 
arguments  did  not  agree  well  with  data.  So  the  actual  forms  used 
for  a  and  S  are  empirically  based.  However,  the  basic  physical 
reasoning  for  the  stability  parameterization  is  as  follows.  In  the 
lower  limit  as  a  approaches  zero,  flow  adjustments  at  each  level 
become  two  dimensional  (as  in  WOCSS) .  This  limit  is  physically 
consistent  with  very  strong  stability.  In  the  upper  limit  as  a 
approaches  unity,  the  solution  becomes  consistent  with  solutions 
for  potential  flow.  So  the  purview  of  NUATMOS  seems  to  be  for 
neutral  to  stable  cases.  Indeed,  for  a  =  1  NUATMOS  was  found  to 
agree  closely  with  exact,  analytic  potential  flow  solutions  which 
are  available  for  idealized  perturbations  such  as  hemispheres, 
half-cylinders,  ellipsoids,  and  polynomial  hills  imbedded  within  an 
otherwise  featureless  horizontal  plane  (Ross  et  al.  1988).  Such 
solutions  do  not  include  adjustments  due  to  thermal  effects  or 
turbulent  drag.  NUATMOS  also  does  not  include  upslope  flow 
adjustments,  ad  hoc  or  otherwise. 
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Dispersion  from  the  puff  model  CITPUFF  with  the  flow  field  supplied 
by  NUATMOS  was  also  tested  against  data  from  the  Cinder  Cone  Butte 
study  (Ross  et  al.  1991)  .  Cinder  Cone  Butte  is  a  lone  roughly 
axisymmetric  hill  rising  from  an  otherwise  nearly  flat  plain. 

The  available  documentation  is  not  clear  on  this  issue.  However, 
in  principle  NUATMOS  seems  to  demand  an  initially  uniform 
undisturbed  flow  field  upon  which  to  make  its  stability  and  terrain 
based  adjustments.  Isolated  hills  such  as  Cinder  Cone  Butte  are 
ideally  suited  to  meet  such  a  demand.  However,  obtaining  a 
pristine  upstream  vector  may  be  difficult  in  Vandenberg's  highly 
complex  terrain.  A  tower  averaged  wind  vector  will  include  myriad 
terrain  effects.  Overwater  vectors  from  NOAA  buoys  do  not  include 
the  effects  of  surface  roughness,  stability,  or  boundary  layer 
height  changes  seen  over  neighboring  but  otherwise  uniform  land 
surfaces.  The  only  tower  in  a  low  region,  009,  is  subject  to 
Lompoc  Valley  winds  which  are  often  not  representative  of  the 
remainder  of  the  flow  over  south  Vandenberg.  An  objectively 
analyzed  windfield  supplied  from  tower  and  SODAR  data  would  already 
contain  flow  adjustments  which  would  be  repeated  in  NUATMOS.  The 
net  result  might  be  exaggeration  or  distortion  of  certain  flow 
features.  Tests  against  Vandenberg  tower  winds  would  be  required 
to  make  a  proper  evaluation  of  NUATMOS.  We  cannot  suggest  that  the 
Cinder  Cone  Butte  data  are  sufficient  for  validation  purposes. 

CITPUFF: 

The  following  analysis  is  based  on  a  1985  publication  which  may  be 
somewhat  out  of  date.  An  updated  analysis  is  in  order,  if  more 
recent  references  become  available.  CITPUFF  appears  to  be  a 
standard  puff  model  which  includes  three  options  for  calculating 
puff  dispersion;  PG  stability  classes  or  Briggs  functions  for  ct  j, 
as  functions  of  stability  class  and  downwind  distance,  or  a  simple 
22.5°  arc  assumption.  Briggs-type  formulas  are  used  for  plume  rise 
and  effective  stack  height.  The  available  CITPUFF  documentation 
makes  no  mention  of  the  many  add-on  features  demanded  of 
operational  diffusion  models,  such  as:  partial  inversion 
penetration,  puff-splitting,  near-field  deformation,  building  down- 
wash,  vertical  shear  effects,  wet/dry  deposition,  multiple  sources, 
source  strength  calculations,  dense  gas  treatment,  chemistry, 
graphical  user  interfaces,  graphical  output,  computer  hardware 
requirements,  or  CPU  run-time  information  on  various  systems. 
Thus,  many  features  including  MARSS  compatibility  and  current  state 
of  support  for  this  model  are  unknown.  Some  speculation  on  the 
effect  of  terrain  induced  streamline  distortion  on  puff  diffusion 
is  given,  but  has  no  immediate  bearing  on  CITPUFF,  which  uses  flat 
terrain  assumptions. 

CITPUFF  was  compared  with  the  COMPLEX  I  and  COMPLEX  II  models  for 
idealiz,ed  flat  terrain  and  isolated  hill  cases.  Since  in  this 
comparison  CITPUFF  used  a  flow  model  other  than  NUATMOS,  little  of 
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relevance  can  be  inferred  from  this  study.  However,  the  Cinder 
Cone  Butte  study  compared  NUATMOS/CITPUFF  against  models  such  as 
NEUTRAL  and  LIFT  for  flow  above  the  CSH,  and  IMPINGEMENT  and  WRAP 
for  flow  below.  LIFT  and  WRAP  (Venkatram,  1988)  were  developed  as 
improvements  to  the  NEUTRAL  and  IMPINGEMENT  models,  respectively. 
In  turn  these  are  variants  of  the  well  known  CTDM  and  CTDM+ 
modeling  systems  (Strimaitis  et  al.  1983,  1987).  With  respect  to 
scatterplots  of  predicted  versus  observed  concentrations,  NUATMOS/ 
CITPUFF  appeared  to  be  at  least  as  accurate  as  the  best  of  the 
other  models,  namely,  NEUTRAL  and  WRAP.  In  fact  NUATMOS/CITPUFF 
showed  no  statistically  significant  fractional  bias  and  a 
normalized  mean  square  error  significantly  smaller  than  WRAP  in 
these  tests.  Inclusion  of  the  a  adjustment  factor  in  NUATMOS 
caused  a  large  improvement  in  results.  The  authors  also  note  that 
NEUTRAL  performed  much  better  than  its  supposed  improvement,  LIFT. 
However,  as  in  the  NUATMOS  review,  we  do  not  recommend  applying 
results  from  an  isolated  hill  like  Cinder  Cone  Butte  to  Vandenberg, 
one  problem  being  specification  of  a  true  background  mean  wind  in 
the  context  of  complex  terrain.  NUATMOS  seems  to  be  particularly 
sensitive  to  this  issue.  Results  from  a  forthcoming  Tracy  power 
plant  study  using  NUATMOS/CITPUFF  may  be  more  relevant. 

SUMMARY: 

In  summary  NUATMOS  provides  a  simple  way  to  include  the  effects  of 
stability  and  dynamic  wind  speed-up  in  an  otherwise  purely  mass 
conserving  model.  Yet,  the  algorithm  is  empirical  because  attempts 
to  base  the  a  tuning  parameter  on  physical  considerations  were  not 
successful.  Being  a  recent  development,  NUATMOS  does  not  appear  to 
be  fully  mature.  Validation  studies  in  truly  complex  terrain  are 
not  yet  available  in  the  literature  and  NUATMOS  still  lacks  many 
niceties  found  in  operational  versions  of  diagnostic  models. 
CITPUFF  appears  to  be  a  standard  puff  model  with  all  the  typical 
issues  including  the  puff /plume  averaging  time  and  spectral 
filtering  problem  mentioned  a  number  of  times  above.  If  CITPUFF 
also  has  not  been  updated  yet  to  include  a  number  of  desired 
operational  features,  we  cannot  recommend  further  evaluation  of 
this  model  until  more  relevant  validation  studies  and  better 
documentation  are  available.  This  may  occur  too  late  for  the 
Vandenberg  modeling  itinerary. 
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WOCSS/Adaptive  Plume  (F.L.  Ludwig,  et  al.,  SRI  Int.) 


INTRO: 

WOCSS  (winds  on  critical  streamline  surfaces)  is  a  windflow  model 
from  SRI  International  which  employs  critical  streamline  heights 
and  mass  conservation  but  no  dynamics,  i.e.,  no  momentum 
conservation.  An  initial  guess  windfield  is  supplied  by  available 
combinations  of  tower,  SODAR,  rawinsonde,  and  National  Weather 
Service  based  wind  data.  These  are  1/R"  interpolated  to  the  grid, 
n  is  a  user-specifiable  variable,  usually  between  one  and  three. 
The  coordinate  system  follows  flow  surfaces  defined  curvilinear ly 
by  critical  streamlines.  Flow  is  partitioned  into  several  non¬ 
interacting  layers,  each  of  which  are  assumed  to  have  constant 
temperature  lapse  rates.  Because  the  surfaces  are  decoupled, 
recirculations  from  slope  winds,  valley  winds,  and  seabreezes  which 
depend  on  vertical  flow  cannot  be  modeled.  This  also  applies  to 
the  typical  secondary  circulation  around  stratus  cloud  edges  as 
discussed  in  the  CALMET  review. 

PHYSICS: 

The  local  height  of  each  surface  is  defined  by. 
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where  Zjq  is  the  height  of  the  bottom  of  the  jth  flow  surface  taken 
over  the  lowest  terrain  in  the  domain.  H,  and  are  the 
local  terrain  height  and  the  highest  and  lowest  terrain  heights  in 
the  domain,  respectively;  f  is  an  adjustable  function.  Using  the 
standard  Froude  number  concept,  the  maximum  height  of  the  flow 
surface  in  the  domain  is  given  by 


^y.max  ~  ^j.O  ^J.O 


g_ 

^o'^O 


-112 


(4) 


Vjo  is  the  windspeed  for  level  j  at  the  lowest  point  in  the  domain. 
In  some  locations  we  may  have  the  condition,  z  <  H,  so  the 
coordinate  system  defined  by  the  flow  surfaces  will  then  intersect 
the  terrain.  With  divergence  minimization,  this  forces  the  wind  to 
flow  around  some  terrain  obstacles  under  stable  conditions  because 
it  lacks  the  kinetic  energy  to  surmount  obstacles  taller  than 
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As  with  the  other  CSH  based  models,  the  underlying  Froude  criterion 
is  physically  suspect,  since  contrary  to  data  it  implies  zero 
windspeed  at  the  CSH  and  local  minimums  atop  hills.  Smith  (1990) 
showed  that  pressure,  not  potential  energy  variations,  control  wind 
speeds.  Yet,  perhaps  fortuitously,  applications  of  the  streamline 
concept  are  still  in  reasonable  agreement  with  data  (see  CALMET 
review  for  more  detail) . 

Also,  specifying  Vjq  may  be  a  problem  for  Vandenberg,  since  the 
lowest  lying  towers  are  in  the  Santa  Ynez  River  Valley.  The  valley 
towers,  in  particular,  tower  009,  are  subject  to  channeled  down- 
valley  drainage  winds  under  stable  conditions  which  may  be  quite 
different  from  the  domain  scale  flow. 

Unlike  the  primitive  3-D  continuity  equation  used  in  CALMET,  WOCSS 
accelerates  the  computation  by  integrating  the  vertical  divergence 
across  finite  layers  whose  local  thicknesses,  $j,  are  defined  as 
fixed  fractions  of  the  local  CSH.  That  is,  the  primitive  3-D 
continuity  equation  can  be  re-arranged  as, 

dw  _  .  du  ^  dv , 
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From  (5)  we  obtain  the  vertical  velocity  change.  Aw,  across  a  layer 
by  parsing  the  vertical  divergence  and  layer  integrating  to  give. 


Aw  =  w,  -  w, 
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Note  that  eqn.  (6)  uses  layer  averaged  velocity  gradients.  We  can 
also  obtain  the  same  vertical  velocity  change  as  in  eqn.  (6)  by 
subtracting  a  vertical  velocity  equation  applied  at  the  bottom  of 
the  layer  from  a  similai  one  applied  at  the  top. 
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If  we  equate  the  two  expressions  for  vertical  velocity  change  and 
re-arrange,  we  obtain, 
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Employing  the  inverse  of  the  chain  rule  for  differentiation  gives, 
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This  reduces  to 
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if  the  layer  thickness,  *j,  does  not  change  with  time.  So  a  fixed 
#j  replaces  iterative  computation  of  the  vertical  divergence,  aw/az, 
resulting  in  a  considerable  savings  in  computer  time. 

As  a  further  approximation,  WOCSS  defines  the  variables. 


u  *  =  uAz  ,  and  v  *  -  vAz 


(11) 


where,  instead  of  ♦j  Az,  the  average  separation  distance  between 
adjacent  surfaces  is  used.  The  assumption  is  that,  if  the  slopes 
of  the  flow  surfaces  are  not  too  large,  the  continuity  constraint 
can  then  be  approximated  by. 
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dx  ay 


Even  so,  the  divergence  minimization  procedure  is  still  the  time- 
bottleneck.  But  with  the  approximations  used,  tests  indicate  that 
on  new  workstations  WOCSS  can  probably  meet  the  emergency  response 
time,  resolution,  and  domain  size  requirements  stipulated  for 
Vandenberg. 

In  WOCSS  continuity  should  accelerate  winds  over  high  terrain  where 
the  layers  are  thinner.  This  is  consistent  with  tests  which  show 
speedup  over  or  around  hills  under  stable  conditions.  But,  the  CSH 
concept  upon  which  WOCSS  is  based  is  physically  suspect,  so  it  is 
not  surprising  that  there  are  qualitative  differences.  That  is, 
for  vertically  uniform  stable  stratif icarion,  linear  diagnostic 
iTt'  iels  also  produce  a  speedup  over  hills  which  are  driven  by  stress 
interactions  with  local  pressure  gradients  rather  than  continuity 
of  flow  through  thinning  layers.  Unlike  WOCSS,  such  linear  models 
show  a  secondary  downwind  speed  maximum  which  eventually  overwhelms 
the  hilltop  maximum,  as  stability  increases.  Stagnation  at  the 
upstream  base  of  the  hill  is  also  augmented  (Carruthers  and  Hunt, 
1990)  . 

Though  there  is  no  explicit  moisture  in  WOCSS,  optional  cloud  cover 
data  can  be  used  to  estimate  net  radiation  at  the  surface  and  hence 
PG  stability  class.  WOCSS  also  estimates  turbulence  kinetic  energy 
(tke)  and  Richardson  number.  However,  they  are  not  yet  used  in  the 
Adaptive  plume  computations. 

SUMMARY; 

WOCSS  takes  a  more  sophisticated  approach  to  mass  conservation  than 
MELSAR  and  MATHEW  by  parsing  divergence  minimization  into  discrete 
layers  so  that  a  pseudo-2-D  equation  can  be  used.  This  results  in 
large  computational  savings.  However,  the  parsing  is  done  in  terms 
of  the  critical  streamline  height  concept,  long  in  vogue,  but 
recently  shown  to  be  physically  suspect.  It  also  lacks  the 
momentum  conservation  of  LINCOH  or  empirical  approaches  for  dealing 
with  dynamical  momentum  effects  and  slope  flows  found  in  CALMET. 
On  the  other  hand,  we  have  found  that  it  does  run  at  500  meter 
resolution  for  a  25  x  40  km  grid  on  a  386  PC  within  “four  minutes 
using  four  levels.  So  on  a  faster  machine,  we  expect  that  WOCSS 
will  comply  with  the  suggested  time  constraints,  while  we  doubt 
that  CALMET  can. 


MACHWIND:  (R.E.  Meyers,  U.S.  Army  Atmospheric  Sciences  Lab) 

MACHWIND  is  a  derivative  of  the  WOCSS  diagnostic  windflow  model 
which  has  been  adapted  for  use  at  Vandenberg.  MACHWIND  improves 
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computationally  upon  the  original  WOCSS  version  by  using  a  three- 
nested  grid  method  to  reduce  the  time  required  for  divergence 
minimization,  the  most  time-consuming  step.  The  coarse  and  medium 
grids  are  set  at  1/4  and  1/2  the  fine  grid  resolution, 
respectively.  From  an  initial  objective  analysis,  a  minimally 
divergent  windfield  is  then  computed  on  the  coarse  grid. 
Respecting  the  fixed  winds  represented  by  the  tower  sites,  this 
field  is  interpolated  to  the  medium  grid  and  becomes  the  initial 
field  for  divergence  minimization  on  the  medium  grid.  The  process 
is  repeated  for  the  fine  grid.  The  multi -grid  procedure  reduces 
the  nvimber  of  iterations  required  overall  by  eliminating  the 
coarser  scale  divergences  first. 

Among  the  Vandenberg  adaptations  are  the  use  of  SODAR  and  tower 
data  to  produce  ncn-uniform  inversion  base  heights.  The  MACHWIND 
inversion  height  at  a  given  grid  point  is  based  on  the  terrain 
height  averaged  over  a  ten  kilometer  radius,  plus  a  l/r“ 
interpolation  of  boundary  layer  depths  taken  from  the  four  coastal 
SODARs.  So  the  modeled  inversion  would  be  basically  flat,  if  the 
inversion  base  height  is  derived  from  one  SODAR  measurement  plus  a 
terrain  height  averaged  over  a  very  large  radius.  On  the  other 
hand,  if  the  radius  were  quite  small,  then  the  inversion  base  wculd 
become  essentially  terrain  following.  The  default  ten  kilometer 
terrain  averaging  radius  is  an  optimization  based  on  the  following 
consideration.  If  tower  based  temperatures  match  those  within  the 
subsidence  inversion  (as  indicated  by  rawinsondes)  it  is  likely 
that  the  inversion  intersects  the  terrain  at  such  sites.  It  was 
found  that  a  ten  kilometer  averaging  radius  placed  the  maximum 
number  of  such  towers  within  the  inversion.  We  find  this  procedure 
somewhat  ad  hoc  but  perhaps  serviceable.  Kamada  et  al.  1989  and 
1990,  showed  for  the  summer  day-time  period  that  the  highest 
correlations  with  mobile  SODAR  data  at  Vandenberg  were  obtained 
assuming  a  terrain  parallel  inversion  height.  Semi-parallel 
inversion  height  assumptions  showed  slightly  lower  correlations. 

Already  planned  for  future  versions  of  MACHWIND  is  an  estimate  of 
the  local  temperature  lapse  rates  to  improve  both  the  windspeed 
interpolation  between  vertical  levels  (now  assumed  logarithmic)  and 
to  also  allow  the  sigma  surfaces  to  be  based  on  hydrostatic  iso¬ 
pressure  surfaces,  rather  than  terrain.  We  agree  with  Hanna  in 
Venkatram  and  Wyngaard,  1988,  that  better  lapse  rates  are  likely  to 
be  obtained  from  surface  energy  balance  methods  than  tower  data. 
The  thermister  accuracy  required  for  the  tower  method  requires 
calibrations  too  frequent  and  detailed  for  routine  operational  use. 
A  sample  surface  energy  balance  method  is  given  in  Appendix  A  of 
the  LINCOM/RIMPUFF  Mt.  Iron  Comparison  (Kamada,  1992) . 

To  accelerate  the  incorporation  of  tower  data  in  the  divergence 
minimization  process,  a  time  transient  velocity  equation  has  been 
suggested  which  includes  advection  and  a  pseudo-pressure  gradient 
force  based  on  the  divergence.  One  suggested  form  of  the  pressure 
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gradient  force  is  based  on  the  deep  convection  form  of  the 
continuity  equation  which  includes  density  changes.  Even  with  time 
tendencies  of  density,  sound  waves  will  not  arise,  since  the  time 
tendency  of  vertical  velocity  is  not  included.  The  deep  convection 
continuity  equation  is  used  in  prognostic  models  such  as  RAMS  to 
simulate  cumulus  circulations  which  involve  large  fractions  of  the 
troposphere.  The  second  suggested  form  for  the  pseudo-pressure 
gradient  force  involves  a  large  coefficient  times  the  divergence  of 
the  velocity  divergence.  These  suggested  forms  seem  rather 
involved,  considering  that  the  purpose  seems  to  be  to  simply 
diffuse  data  induced  divergence  efficiently.  Moreover,  the  use  of 
time  dependent  equations  demands  considerably  more  computer  power. 
We  question  whether  such  a  set  can  meet  Vandenberg's  emergency 
response  constraints  on  currently  available  workstations. 


SUMMARY : 

MACHWIND  is  a  derivative  of  the  above  reviewed  WOCSS  model  which 
has  been  adapted  to  use  at  Vandenberg.  The  main  advantages  over 
WOCSS  seem  to  be  the  faster  divergence  minimization  scheme  and  non- 
uniform  inversion  base  heights.  .Also  planned  are  a  conversion  to 
hydrostatic  sigma  surfaces  and  the  use  of  time  transient  truncated 
Navier-Stokes  equations  to  facilitate  the  diffusion  of  data  induced 
divergence  into  the  windfield.  The  latter  would  make  MACHWIND  a 
considerably  different  model  from  WOCSS.  However,  we  question  the 
latter  implementation  for  Vandenberg  emergency  response,  given 
currently  available  workstation  speeds.  With  appropriate  lapse 
rate  estimates,  it  appears  that  MACHWIND  ruay  be  the  most  suitable 
of  the  current  models  designed  for  stable  flow  conditions,  as  well 
as  the  one  nearest  to  operational  status  at  Vandenberg. 

Adaptive  Plume  Model: 

PHYSICS: 

Rather  than  a  spherical  puff,  the  adaptive  plume  model  describes  a 
Lagrangian  volume  using  five  different  points  in  the  vertical  with 
an  initial  separation  of  three  meters.  All  material  within  the 
puff  lies  between  the  highest  and  lowest  points.  The  fraction  of 
total  material  between  any  two  of  the  five  points  is  assumed  to  be 
given  by  a  piecewise  Gaussian  function.  The  resulting  nearly  five¬ 
fold  increase  in  computational  demands  is  somewhat  offset  by  an 
attempt  to  merge  puffs  within  a  distance  of  ay/2  as  well  as  purge 
puffs  which  have  exited  the  domain.  The  purpose  of  the  five  points 
is  to  allow  for  vertical  differences  in  advection  and  diffusion 
presented  by  the  windfield  model  which  may  be  more  complex  than 
that  given  by  standard  similarity  profiles.  Thus,  vertical  shear, 
plume  rise,  and  the  effect  of  growing  internal  boundary  layers  all 
may  be  treated  within  the  model.  The  very  small  initial  separation 
between  points  suggests  that  the  near-field  slug/puff  adaptation 
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used  by  CALPUFF  is  not  necessary  for  the  Adaptive  plume  model. 
However,  the  need  to  release  small  puffs  frequently  does  increase 
computational  cost. 

Plume  rise  is  handled  via  adjustments  of  the  standard  Briggs 
equations  applied  separately  to  each  of  the  five  points.  100 
receptor  sites  may  be  specified.  At  present  vertical  dispersion  is 
based  on  PG  stability  classes,  which  we  noted  above  tend  to 
underestimate  near-field  vertical  diffusion  for  surface  releases 
under  convective  conditions.  Complete  specular  reflection  is 
assumed  at  the  surface,  such  that  deposition  effects  are  not 
computed.  Puffs  are  not  allowed  to  split  and  there  is  no  provision 
for  treating  mviltiple  sources  within  a  single  run.  As  with  CALPUFF 
there  do  not  appear  to  be  provisions  for  spectral  filtering  of  the 
dispersion  coefficients.  So  adjustments  would  be  required  when 
using  this  model  with  a  prognostic  or  frequently  updated  diagnostic 
wind  field.  Graphical  output  is  now  available  in  one  version  of 
the  model  and  it  has  been  run  on  MARSS  compatible  DEC  systems. 

VALIDATION; 

The  model  was  tested  against  LIDAR  data  collected  on  May  5  and  15, 
1980  from  the  EPRI  project  at  the  Kincaid  power  plant  in  Illinois. 
Good  topological  agreement  was  obtained  between  the  model  and 
vertically  oriented  2-D  frames  from  the  LIDAR.  Thus  far, 
comparisons  of  surface  footprint  data  from  other  studies  are  not 
available. 


WOCSS/Adaptive  Plume  SUMMARY: 

The  adaptive  plume  concept  appears  to  have  some  merit  and  invites 
further  study.  However,  the  model  is  in  research  level,  proof -of- 
concept  status  rather  than  in  a  fully  operational  system  mode  which 
includes  all  features  necessary  for  dispersion  calculations.  Like 
the  other  reviewed  models,  relevant  chemistry,  source  terms  for 
liquid  spills,  and  options  for  treating  dense  gases  are  not 
available.  Near-field  building  down-wash,  multiple  sources, 
overwater  transport,  and  wet/dry  deposition  effects  are  also  not 
implemented. 

Both  WOCSS  and  the  Adaptive  plume  models  are  in  the  public  domain, 
although  some  restrictions  apply.  The  code  is  modular  and 
modifiable.  Technical  support  is  available  from  SRI,  but  is  not 
likely  to  come  gratis.  Versions  exist  in  VMS  Fortran,  so  MARSS 
compatibility  is  probably  not  an  issue.  In  fact,  most  of  the 
reviewed  models  have  been  written  in  ANSI  F77  standard  Fortran. 
Thus,  portability  should  be  a  trivial  matter  in  most  cases.  A 
graphical  output  format  is  being  developed  at  San  Jose  State  ,  but 
there  are  no  current  plans  for  a  full  graphical  user  interface 
(F.L.  Ludwig,  personal  communication,  1992). 
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For  satisfactory  operation  at  Vandenberg,  the  Adaptive  plume  model 
would  probably  require  extensive  additions  to  the  basic  algorithms 
to  include  the  many  add-on  features  required  of  operational  models. 
Given  the  right  personnel,  this  might  be  done  in  less  than  one  man- 
year.  However,  we  suggest  that  other  more  mature  models  be  studied 
before  attempting  such  an  undertaking. 
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RISO  National  Labs, 


LINCOM-RIMPUFF:  (Troen,  Mikkelsen,  et  al., 

Denmark) 

LINCOM  (I.  Troen  et  al.): 

INTRO: 

LINCOM  belongs  to  the  only  class  of  diagnostic  model  which  retains 
both  mass  and  momentum  conservation.  The  first  proponents  were 
Jackson  and  Hunt  (1975),  followed,  e.g. ,  by  Mason  and  Sykes  (1979), 
Walmsley  et  al.  (1982),  Taylor  et  al.  (1983),  Beljaars  et  al. 
(1987),  and  Hunt  et  al .  (1988a).  LINCOM  uses  linearized  versions 
of  the  perturbation  momentum  equation  with  first  order  turbulent 
diffusion,  the  standard  shallow  convection  form  of  mass  continuity, 
and  equation  of  state.  Thus,  to  first  order,  LINCOM  conserves 
momentum  as  well  as  mass  by  treating  the  effects  of  pressure  force 
dynamics,  turbulent  stress,  advection,  coriolis,  and  continuity, 
Currently,  the  thermodynamic  energy  or  temperature  equation  is 
neglected.  So,  buoyancy  induced,  non-neutral  conditions  are  not 
treated  explicitly  in  the  governing  equations.  (A  thermal  version 
of  LINCOM,  LINCOM-T  is  under  development.) 

Actually,  in  version  1.0a  post-process  objective  analysis  anchors 
the  wind  field  to  match  the  tower  vectors  exactly.  However,  in 
this  process  true  mass  conservation  is  lost.  The  objective 
analysis  uses  a  terrain  modified  1/R^  interpolation,  wherein  each 
tower's  circle  of  influence  is  warped  into  ovals  matching  the 
aspect  ratio  of  the  local  terrain  contours  surrounding  the  tower 
site.  This  should  be  an  improvement  over  the  standard  Barnes 
technique.  If  only  surface  data  is  supplied,  LINCOM  extrapolates 
vertically  on  the  basis  of  similarity  theory  to  obtain  a  mean 
background  wind.  Version  1.0b  is  significantly  different  in  that 
it  retains  mass  conservation  by  using  a  fitting  techique  which  does 
not  attempt  an  exact  match  to  tower  wind  vectors.  Moreover,  for  a 
given  domain,  orthogonal  basis  functions  for  the  wind  vectors  are 
pre-computed  and  stored  on  hard  disk.  So  assembly  of  the  final 
wind  field  only  involves  a  non-iterative  variational  search  for  the 
domain  mean  wind  vector,  which  when  added  to  the  pre-computed 
perturbation  field,  provides  the  best  fit  to  the  available  data. 
Thus,  this  model  may  be  one  to  two  orders  of  magnitude  faster  than 
the  other  windfield  models  reviewed  here. 

PHYSICS: 

Both  LINCOM  versions  solve  the  equations  by  Fourier  transforming 
the  terrain  as  well  as  the  flow.  For  each  wave  number  in  the 
spectral  domain,  pre-calculated  analytical  solutions  are  available 
for  linearized  equations.  Unlike  the  divergence  minimization 
schemes,  time-consuming  iterative  numerical  procedures  are  not 
involved.  Analytic  solutions  could  also  be  obtained  for  individual 
grid  points  in  finite  difference  mode.  The  reason  for  the 
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transformation  from  grid  points  to  wave  numbers  is  that  only  "1/3 
as  many  wave  numbers  are  required  as  grid  points  in  each  dimension. 
This  reduces  computational  requirements  by  an  order  of  magnitude 
for  a  2-D  transform  (minus  the  relatively  small  time  involved  in 
computing  the  transform  and  inverse  transform) . 

Though  useful  in  setting  up  a  domain,  in  LINCOM  version  1.0  the 
speed  gained  by  using  Fourier  transforms  is  actually  academic 
because  for  operational  purposes  the  solutions  are  pre-calculated . 
Thus,  version  1.0  stores  72  pre-calculated  solutions  for  the 
terrain  induced  perturbations  for  a  given  domain  on  hard  disk  in 
three  degree  increments  around  360°  of  arc  for  reference  1  ms  '  mean 
background  winds.  Version  1.0  then  retrieves  a  stored  perturbation 
windfield  based  on  the  mean  wind  direction  supplied  by  that  tower. 
Since  flow  above  the  surface  layer  is  assumed  to  be  frictionless, 
the  perturbation  field  vectors  are  simply  multiplied  by  the  tower 
wind  magnitude  and  added  to  the  mean  vector  to  produce  the  final 
field  based  on  that  particular  tower.  Version  1.0  then  combines 
separate  fields  obtained  for  each  tower,  using  the  terrain  modified 
1/R-  interpolation.  Since  the  perturbations  are  constrained  to  be 
zero  at  each  tower  site,  the  final  field  provides  an  exact  match  to 
the  tower  winds,  while  allowing  a  dynamically  based  interpolation 
scheme  between  towers.  This  accounts  for  the  well  validated  wind 
•'speed-up"  in  the  surface  layer  over  hills.  Although  the 
individual  perturbation  fields  are  mass  conservative,  mass 
conservation  is  lost  in  the  combinatory  process.  Version  1.0  was 
used  for  the  Mt.  Iron  tracer  study. 

Another  aspect  of  linear  systems  exploited  in  version  1.0b  is  that 
the  total  flow  is  expressible  as  a  linear  combination  of  orthogonal 
solutions  of  the  perturbation  flow  field.  This  means  that  all  the 
basic  features  of  neutral  flow  can  be  expressed  as  a  linear 
combination  of  just  two  fixed  orthogonal  fields  which  may  be  pre¬ 
calculated  for  any  fixed  set  of  terrain.  So  in  practice  two, 
rather  than  72,  perturbation  fields  are  stored,  based  on  1  ms'  due 
northerly  and  easterly  mean  winds,  a  36-fold  savings.  Instead  of 
combining  multiple  fields  based  on  exact  anchoring  to  the  towers, 
a  single  field  is  obtained  by  variationally  optimizing  the  fit  to 
the  horizontal  tower  winds.  Thus,  this  procedure  yields  both 
momentum  and  mass  conservation  in  a  single  non-iterative  pass.  And 
because  the  essential  flow  features  are  pre-stored  rather  than 
calculated,  it  is  also  very  fast.  On  a  15  MIPS  486/33  PC,  LINCOM/ 
RIMPUFF  required  three  minutes  to  simulate  nine  hours  of  continuous 
plume  dispersion  using  windfields  updated  every  ten  minutes  on  a 
200  x  200  X  5  grid  at  50  meter  resolution.  So  here  LINCOM  requires 
less  than  three  seconds  to  formulate  and  output  each  of  54  wind 
fields . 

In  actual  operation  the  input-to-output  time  for  version  1.0b 
depends  on  the  trivial  fitting  routine  and  hard-disk  retrieval  of 
pre-calculated  results.  Version  l.O  is  generally  somewhat  slower. 
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Both  versions  have  also  been  run  on  a  DEC  VAX  5000  workstation,  so 
software  compatibility  with  MARSS  seems  likely. 

The  complete  linearity  presumed  in  the  pre-calculation  method 
assumes  that  the  flow  is  inviscid  in  the  outer  boundary  layei , 
which  is  less  reasonable  within  the  inversion  (which  LINCOM  assumes 
is  terrain  following) .  LINCOM  also  lacks  time  tendencies  in  the 
velocity,  so  terrain  forced  dynamics  propagate  instantly  over  the 
whole  domain,  rather  than  at  some  realistic  set  of  phase  speeds. 
This  problem  is  shared  by  all  diagnostic  models.  Linearization 
also  means  that  higher  order  perturbation  terms  are  neglected.  In 
a  prognostic  scheme,  such  solutions  would  gradually  diverge  from 
reality.  However,  for  a  diagnostic  model,  where  only  the  present 
time  step  is  considered,  this  is  not  as  significant  an  issue. 
Still,  linear  models  cannot  treat  large  non-linear  synergisms.  One 
example  might  be  the  south  side  of  Honda  Ridge  during  north- 
westerlies  where  roughness  and  stability  changes  occur  along  with 
flow  separation. 

The  vertical  velocity  is  basically  assumed  equal  to  the  horizontal 
velocity  at  that  site  times  the  sine  of  the  terrain  slope.  This 
only  holds  well  for  slopes  of  less  than  twenty  degrees  (Hunt, 
personal  communication,  1990).  However,  ACTA's  analysis  of  the 
Vandenberg  terrain  concludes  that  at  500  meter  resolution,  little 
of  the  terrain  shows  slopes  exceeding  twenty  degrees  (Conley  et 
al.,  1990). 

Transformation  of  a  limited  mesoscale  area  to  the  spectral  domain 
involves  the  assumption  of  periodic  boundary  conditions.  So  wave 
energy  leaving  one  boundary  immediately  re-enters  the  domain  from 
the  opposite  boundary.  Thus,  LINCOM  employs  a  10  km  buffer  zone 
over  which  terrain  is  gradually  relaxed  on  all  sides  to  sea  level. 
This  greatly  damps  wave  motions,  but  not  completely.  Quantitative 
assessment  of  the  influence  of  this  effect  on  flow  is  not  presently 
available . 

Exact  matching  to  tower  winds  is  not  always  desirable,  since  tower 
winds  may  represent  very  local  conditions  which  scarcely  influence 
the  bulk  flow  (see  validation  effort  below) .  This  can  be  regarded 
as  a  form  of  aliasing.  LINCOM  1.0  also  presumes  that  each  tower 
contributes  its  version  of  the  undisturbed  mean  flow  vector.  But 
the  in-terrain  towers  are  all  influenced  by  the  surrounding  complex 
terrain  and  can  hardly  be  considered  representative  of  undisturbed 
background  flow.  For  this  reason  Vandenberg  tower  052,  sited 
within  a  reasonably  flat  mesa  portion  of  the  domain,  was  used  to 
represent  background  flow  in  the  Mt.  Iron  comparison.  LINCOM  1.0b 
places  less  importance  on  the  meaning  of  the  background  wind  by 
deriving  it  through  empirical  fitting.  This  seems  more  honest. 

Version  1.0b  also  does  not  suffer  the  loss  of  mass  conservation  or 
the  local  aliasing  problem.  However,  its  tower  fitting  technique 
cannot  yield  local  slope  flows  or  stability  effects,  since  they  are 
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not  treated  in  the  governing  equations.  Hence,  the  fitting 
technique  in  version  l.Ob  may  produce  poor  results  when  such 
effects  are  important  in  determining  local  wind  directions. 
Stability  alters  local  winds  strongly  during  drainage  flows.  So 
Froude  balance  models  might  be  more  appropriate  than  neutral  LINCOM 
under  strongly  stable  conditions.  Real  slope  flow  forcing 
increases  as  the  square  of  the  slope's  length,  so  small  scale 
slopes  do  not  contribute  much  to  slope  flows.  For  example  the 
Vandenberg  seabreeze  is  often  augmented  significantly  by  a 
southwesterly  up”slope  flow.  This  flow  is  due  to  heating  of  the 
large  scale  slope,  extending  from  the  coast  to  the  Santa  Ynez 
ridge,  about  75  to  100km.  However,  for  LINCOM's  purposes  the 
domain  mean  background  wind  incorporates  both  this  effect  as  well 
as  the  Seabreeze. 

Both  LINCOM  versions  assume  a  terrain  parallel  inversion.  This 
agrees  generally  with  daytime  data  from  the  Vandenberg  area  studies 
of  Kamada  et  al.  (1990a,  b) .  However,  such  an  assumption  may  not 
be  valid  under  strongly  synoptic  or  stable  conditions. 

VALIDATION: 

Version  1.0  has  been  tested  together  with  RIMPUFF  for  eight 
representative  cases  from  the  1966-67  Mt.  Iron  tracer  study 
conducted  by  Battelle  Corp.  over  south  Vandenberg.  In  some  cases 
modeled  plumes  tracked  measured  dosages  quite  well,  due  in  part  to 
"tower  anchoring".  In  other  cases,  the  real  plume  floated  above 
the  local  canyon  flow  represented  by  towers.  A  simple  improvement 
in  version  1.0  would  be  to  alter  the  1/R^  dependence  used  in  the 
objective  analysis  to  1/R",  where  n  depends  on  local  terrain  slope 
and  stability o  This  would  tend  to  limit  the  influence  of  highly 
localized  effects  on  the  overall  wind  field  and  also  limit  the 
degradation  of  mass  conservation.  However,  this  adjustment  has  not 
been  implemented  or  tested  as  yet. 

The  surface-] ayer,  ridge-top  "speed-up"  recognized  in  LINCOM's 
physics  already  tends  to  localize  some  tower  influences.  That  is, 
many  of  the  Vandenberg  towers  are  sited  along  ridge  tops  rather 
than  at  the  bottom  of  or  along  hills.  Since  ridge  tops  represent 
only  a  small  fraction  of  the  terrain,  these  towers  often  record 
atypically  high  winds.  If  a  flow  model  is  unaware  that  "speed-up" 
is  confined  to  hill  tops,  then  its  objective  analysis  step  will 
over-accelerate  the  rest  of  the  flow  to  match  the  ridge  top  data. 
Such  overly  high  speed  estimates  will  lead  to  overly  thin,  overly 
long  plume  estimates.  Unlike  other  models  LINCOM  should  be  able  to 
avoid  this  effect  during  neutral  and  unstable  conditions. 

Since  in  operation  pre-calculated  solutions  are  used,  the  speed 
gained  by  finding  wave  number  rather  than  grid  point  solutions  is 
irrelevant.  So  it  also  seems  that  the  periodic  boundary  condition 
problem  could  be  avoided  and  nesting  within  larger  scale  models 
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could  be  facilitated  by  instead  storing  finite  difference  grid 
point  solutions  with  zero  gradient  lateral  boundary  conditions. 
This  point  is  being  considered  in  developing  a  thermal  version  of 
LINCOM,  LINCOH-T. 

SUMMARY: 

Given  the  500m  resolution  and  50x50  km  domain  size  demonstrated  for 
Vandenberg,  LINCOM  appears  to  be  reasonably  accurate,  v/hile  also 
capable  of  running  in  near  real-time.  LINCOM  can  easily  run  using 
on-line  data  updated  every  minute.  It  is  also  considerably  more 
sophisticated  than  MELSAR  and  less  ad  hoc  than  CALMET  because  it 
belongs  to  a  class  of  linear  models  which  treat  neutral  flow 
physics,  including  the  turbulent  drag  and  pressure  gradients 
neglected  in  pure  mass  conservation/ Froude  number  based  models. 
I.e.,  LINCOM  belongs  to  the  only  class  of  diagnostic  models  which 
includes  both  mass  and  momentum  conservation.  Since  many 
Vandenberg  towers  are  sited  at  ridge  tops,  a  major  practical 
advantage  of  such  models  is  cognizance  of  ‘'speed-up”  over  hills. 
In  reality  this  "speed-up"  occurs  not  only  for  neutral  but  also 
extends  to  both  unstable  and  stable  conditions.  So  in  its 
objective  analysis,  LINCOM  should  avoid  over-estimates  of  speed  in 
lower  terrain,  i.e.,  the  bulk  of  the  domain. 

Among  the  class  of  linear  dynamics  models  the  chief  advantage  cf 
LINCOM  is  the  pre-calculation  of  purely  analytic  solutions.  This 
makes  it  many  times  faster  than  other  models  of  this  type  (e.g., 
FLOWSTAR,  1988) .  The  chief  disadvantage  is  that  any  non-neutral 
physics  is  currently  treated,  not  in  the  governing  equations,  but 
by  objectively  analyzed  adjustment. 

The  basis  for  current  work  on  extending  LINCOM  to  LINCOM-T  is  that 
in  principle  it  is  possible  to  treat  stable  and  unstable  cases 
v/ithin  the  linear  context  analytically.  Hunt  maintains  that  this 
may  be  true  provided  that  the.  hill  Froude  number,  F^  s  UqH/N  >  1. 
Here  Ug  is  mean  wind  speed,  H  is  hill  height  from  1/2  the  hill 
width,  and  N  =  [  (g/Q)  (5e/3z)  ] is  Brunt-Vaisala  frequency.  For  Fy 
<  1,  the  Hunt  group  claims  that  air  will  tend  to  be  blocked  by  or 
flow  around  a  hill  in  which  case  Froude  balance  models  are  probably 
more  appropriate  (Carruthers  and  Hunt,  1990) .  However,  as 
discussed  in  th.^  CALMET  review.  Smith's  recent  work  suggests  that 
the  kinetic/potential  energy  balance  which  forms  the  basis  for  the 
standard  Froude  number  criterion  is  false  (Smith,  1990) .  If  so, 
LINCOM-T  may  apply  eventually  to  all  stability  conditions. 

Until  LINCOM-T  is  developed,  models  such  as  MACHWIND  may  be  the 
best  currently  available  for  such  stable  conditions.  So  for  now, 
these  two  types  of  models  could  be  woven  easily  into  a  selection 
system  based  on  measured  or  estimated  atmospheric  stability. 
However,  at  Vandenberg  strongly  stable  flows  are  not  nearly  as 
common  as  moderately  stable,  near-neutral,  and  unstable  ones,  due 
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to  the  influence  of  the  marine  boundary  layer, 
daytime  when  cold  spill  potential  is  higher  due 
exception  is  larger  scale  drainage  out  of  the 
In  winter  such  down-valley  drainage  may  not 
morning  in  the  face  of  a  weak  seabreeze. 


particularly  during 
to  operations.  One 
Santa  Ynez  valley, 
reverse  until  mid- 
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RIMPUFF  (T.  Mikkelsen  et  al.): 
PHYSICS: 


RIMPUFF  employs  the  LINCOM  flow  field  to  advect  contaminant  puffs 
downstream.  Puff  growth  is  controlled  by  local  turbulence  level, 
using  spectrally  filtered  relative  diffusion  in  time  and  space.  As 
discussed  in  the  analysis  of  CALMET,  the  filtering  is  undertaken  to 
eliminate  eddy  scales  which  contribute  principally  to  puff  meander 
rather  than  diffusion.  However,  if  the  mean  flow  is  not  updated 
frequently,  RIMPUFF  re-introduces  the  larger  eddy  scales  to 
simulate  stochastic  meander  of  the  puff  centers  of  mass.  This  is 
particularly  important  when  puff  progeny  are  closely  spaced  shortly 
after  puff  splitting.  RIMPUFF  initializes  100m  puffs  and  allows 
them  to  split  horizontally  when  their  diameter  exceeds  the  500m 
LINCOM  grid  spacing  or  vertically  under  shear  conditions.  The 
practical  limit  on  PC  based  computers  is  several  hundred  puff 
progeny.  Lateral  diffusion  depends  on  tower  based  parametrizations 
of  the  lateral  dispersion  coefficients,  These  are  a  function 
of  averaging  time,  stability,  and  wind  direction.  For  vertical 
diffusion,  RIMPUFF  currently  relies  on  PG  stability  classes 
modified  for  complex  terrain.  Note  our  above  caveat  on  the  use  of 
PG  classes  for  estimating  plume  behavior  in  the  vertical  under 
convective  conditions. 

Standard  Briggs-type  equations  are  used  in  the  plume  rise  module. 
100%  reflection  is  assumed  at  a  user  specified  inversion  height 
which  parallels  the  terrain.  Dry  deposition  from  the  standard 
source  depletion  method  and  wet  deposition  with  washout  is 
included.  Like  the  other  reviewed  models,  relevant  chemistry, 
source  terms  for  liquid  spills,  and  options  for  treating  dense 
gases  are  not  available.  Near-field  building  down-wash  is  also  not 
implemented.  Graphical  along  with  tabular  output  is  provided,  and 
a  full  graphical  user  interface  is  currently  available  for  UNIX 
based  workstations.  The  source  code,  written  in  ANSI  F77  standard 
Fortran  is  available  to  Vandenberg,  modular,  and  modifiable. 
RIMPUFF  has  been  run  on  DEC  computers  as  well  as  PCs  and  should 
thus  be  MARSS  compatible.  Typical  run-times  are  on  the  order  of 
one  minute  on  a  486/33  PC. 

VALIDATION: 

Comparisons  of  LINCOM/RIMPUFF  with  eight  representative  cases  from 
the  Mt.  Iron  zinc  sulfide  tracer  study  (Kamada  et  al.  1992a,  b) 
indicate  fairly  good  agreement.  However,  modeled  plumes  were 
longer  than  measured  in  seven  of  the  eight  cases.  This  agrees  with 
other  modeling  comparisons  for  Mt.  Iron  data  (Kunkel  and  Izumi, 
1990  and  Kunkel,  1991)  and  results  from  the  Vandenberg  Lompoc 
Valley  Diffusion  Experiment  using  inert  SF^  tracer  in  1989 
(Skupniewicz  et  al.  1991a, c,  1992).  Briggs  (1988)  suggests  that, 
before  the  deployment  of  field  GC  detection  equipment  for  SF^ 
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around  1970,  all  previously  used  tracer  compounds  had  significant 
deposition  velocities.  Carruthers  and  Hunt  (1990)  suggest  that 
surface  concentrations  should  be  independent  of  deposition  velocity 
within  a  downwind  range  given  by  5  z^U/u, .  This  is  just  a  few 
hundred  meters  for  Mt.  Iron,  since  tracer  was  released  from  a 
height,  z^,  of  only  twelve  feet.  So  the  range  discrepancy  seems  due 
to  a  field  study  rather  than  modeling  artifact. 

Case  55  from  Mt.  Iron  shows  immediate  lofting  of  the  plume  under 
convective  conditions,  wherein  RIMPUFF's  dependence  on  PG  classes 
for  vertical  diffusion  may  be  demonstrably  in  error.  Unlike  the 
Diablo  Canyon  study,  the  Mt.  Iron  isopleths  were  relatively  short 
range,  hardly  extending  beyond  five  kilometers.  So  better 
agreement  is  expected  so  long  as  the  initial  wind  vector  measured 
at  the  release  site  is  itself  accurate.  Case  28  of  the  Mt.  Iron 
study  saows  come  discrepancy  in  this  regard.  Much  of  the  disparity 
between  the  model  and  Mt.  Iron  data  occurred  in  the  very  near 
field,  due  to  the  initial  100  meter  puff  radii.  Much  smaller  puffs 
could  have  been  used.  However,  small  puffs  must  be  released 
frequently  to  maintain  accurate  instantaneous  concentrations,  since 
the  distance  between  successive  puff  releases  must  not  exceed  “  1.5 
Uy.  This  constraint  can  be  relaxed  substantially,  if  instantaneous 
concentrations  are  integrated  to  give  dose  exposures,  provided  the 
wind  direction  remains  constant.  In  practice  smaller  initial  puff 
radii  often  compel  increased  computer  time.  This  trade-off  applies 
generally  to  all  puff  models. 

SUMMARY : 

In  summary  RIMPUFF  seems  to  be  a  fairly  sophisticated  model  whose 
chief  advantage  is  that  it  incorporates  spectral  filtering  and 
averaging  time  dependent  dispersion  coefficients  which  are  already 
tailored  for  Vandenberg.  Thus,  RIMPUFF  can  be  run  with  either  a 
static  or  frequently  updated  diagnostic  wind  field,  or  even  a 
prognostic  one.  We  suspect  that  additions  to  RIMPUFF  to  include 
the  AFTOX  and/or  REED-M  source  terms  and  chemistry,  as  well  as  such 
items  as  building  down-wash,  and  similarity  based  vertical 
diffusion  coefficients  would  not  be  as  difficult  as  modifying 
CALPUFF  or  MESOI  for  spectral  filtering.  RISO  Laboratory  does 
support  HEAVYPUFF,  a  1-D  similarity  based  dense  puff  also  driven  by 
LINCOM  (see  following  review) .  However,  HEAVYPUFF  has  not  as  yet 
been  coupled  to  RIMPUFF.  So  with  respect  to  Vandenberg 's  needs  the 
chief  current  disadvantages  of  RIMPUFF  seem  to  be  inaccuracy  in  the 
near-field  due  to  PG  based  vertical  diffusion,  and  the  omission  of 
algorithms  for  estimating  source  terms  and  downwind  chemistry, 
building  dov/n-wash,  cloud/ciear  interface,  and  dense  gas  effects. 
The  LINCOM/RIMPUFF  tandem  is  generally  much  faster  than  the  other 
modeling  systems.  Nine  hours  of  dispersion  involving  wind  fields 
spaced  every  ten  minutes  (54  in  total)  required  only  three  minutes 
on  a  15  MIPS  486/33  PC  for  a  200  x  200  x5  grid  at  500  meter 
resolution. 
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GENERAL  DISCUSSION  OF  DENSE  GASES; 


We  note  that  dense  gas  behavior  may  differ  from  that  of  trace  or 
neutrally  buoyant  gases  in  a  number  of  ways.  For  example,  like 
water  in  air,  the  cloud's  greater  density  will  cause  it  to  slump 
and  spread  laterally  to  produce  a  lower,  broader  cloud,  as  well  as 
induce  it  to  flow  down  hill.  Greater  entrainment  will  occur  at  the 
leading  edge  for  down-slope  flow  due  to  the  enhanced  shear,  if  the 
cloud's  speed  exceeds  ambient  by  a  value  greater  than  the  friction 
velocity. 

If  the  cloud  is  not  cold  relative  to  the  surface  over  which  it 
travels,  its  initial  speed  may  only  slowly  shift  to  ambient  values. 
That  is,  the  large  cloud/air  density  difference  will  suppress  cloud 
top  entrainment  of  ambient  air.  The  strongly  stratified  in-cloud 
density  gradient  will  also  suppress  in-cloud  turbulence  and  thus 
any  side  entrainment. 

On  the  other  hand  the  material  in  question  may  be  stored 
cryogenically  or  released  upon  rupture  of  a  pressurized  tank. 
Sudden  expansion  upon  loss  of  tank  pressure  will  cool  the  cloud.  If 
the  tank  pressure  used  was  sufficient  to  liquify  the  material, 
flashing,  i.e.,  intense  boiling  and  expansive  release  of  a  rapidly 
evaporating  mist  will  occur,  at  the  lower  temperatures  appropriate 
to  ambient  pressure.  In  any  case  the  resulting  cloud  will  be  much 
colder  than  the  underlying  surface.  If  so,  the  surface  will  warm 
the  cloud  convectively ,  thus  opposing  the  turbulence/entrainment 
suppression  induced  by  density.  Material  dilution  during  flashing 
v'ill  also  diminish  any  initial  cloud/air  speed  differences. 

If  the  release  occurs  in  the  wake  of  a  building,  rapid  entrainment 
may  terminate  the  slumping  phase  quite  quickly.  Then  the  cloud 
will  grow  much  more  rapidly  until  it  is  downwind  of  the  wake  or  the 
mean  cloud  height  is  higher  than  surrounding  buildings.  In  any 
event  the  cloud  will  eventually  behave  as  a  neutral  gas  due  to 
sufficient  dilution  with  ambient  air. 

Dense  gas  models  seem  to  fall  into  three  major  categories  of 
complexity:  a)  1-D  models  which  assume  evenly  distributed 
properties  within  a  cylindrical  box;  b)  intermediate  models  which 
assume  similarity  based  vertical  profiles  of  concentration  and 
other  properties  in  the  cross-wind  plane;  and  c)  fully  3-D  Monte 
Carlo  models  which  diffuse  many  particle  packets  of  concentration 
and  other  properties  according  to  mean  advection  velocities 
supplied  by  the  windflow  model,  plus  a  turbulence  simulating, 
random  component. 

The  simple  box  model,  HEAVYPUFF,  is  associated  with  RIMPUFF,  while 
ADPIC  is  representative  of  the  more  complex  Monte  Carlo  approach. 
Both  will  be  reviewed  in  light  of  the  above  differences  between 
neutral  and  dense  gas  behavior. 
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HEAVYPUFF  (N.O.  Jensen,  M.  Nielsen,  RISO  National  Labs,  Denmark) 


HEAVYPUFF  is  a  simple  PC  based,  dense  puff  model  v;hich  is  menu 
driven,  interactive,  and  delivers  graphical  output.  It  requires 
release  inputs  such  as  spill  size,  storage  temperature,  wind  speed, 
ambient  air  and  surface  temperatures,  surface  roughness,  ambient 
Obukhov  length,  L^,  and  whether  the  material  is  flashing  from  liquid 
phase.  Initial  puff  aspect  ratio,  temperature,  and  air/gas  mass 
ratio  are  also  input.  However,  if  flashing  occurs,  the  temperature 
and  mass  ratio  are  coupled  by  enthalpy  and  cannot  be  specified 
independently,  since  the  model  computes  the  amount  of  air  necessary 
to  evaporate  exactly  all  the  material.  Effectively  complete 
evaporation  may  occur  if  depressurization  is  rapid,  or  if  the  heat 
lost  by  flashing  freezes  the  surface  of  the  remaining  liquid  and 
greatly  reduces  the  subsequent  release  rate.  Such  a  release  will 
appear  as  a  sudden  puff  followed  by  a  greatly  diminished  plume. 
But  the  user  must  then  estimate  the  actual  effective  spill  size 
outside  of  HEAVYPUFF. 

Left  unspecified,  the  release  temperature  will  default  to  the 
liquid's  boiling  point  and  the  cloud  height/radius  aspect  ratio, 
h/R,  will  default  to  0.5.  The  integration  time  step  is  user 
specifiable  or  variably  automatic.  Outputs  range  over  a  list  of  20 
variables,  including  centerline  concentrations  and  the  area  in 
which  concentrations  exceed  user  specified  limits. 

The  physics  in  HEAVYPUFF  treats  heat  convection  from  a  warm  surface 
to  a  cold  cloud,  using  surface  layer  similarity  theory.  Ideally, 
this  requires  flat  terrain  and  a  cloud  height  much  greater  than  the 
surface  roughness,  but  much  less  than  the  atmospheric  boundary 
layer  inversion  lAeight.  The  puff  is  modeled  as  a  cylindrical  box 
with  all  properties  evenly  distributed.  Indeed,  over  a  much  warmer 
surface,  a  cold  cloud,  once  convective,  may  tend  to  even 
distributions  due  to  vigorous  turbulent  mixing  induced  by  the 
strong  upward  heat  flux. 

In  HEAVYPUFF  cloud  volume  and  mass  increase  only  by  cloud-top 
entrainment.  A  parametrized,  vertically  averaged  version  of  the 
turbulent  kinetic  energy  (tke)  equation,  similar  to  that  used  for 
atmospheric  boundary  layer  entrainment,  is  used  to  obtain  the 
cloud-top  entrainment  rate.  The  heat  from  the  warmer  air  and  shear 
generated  tke  is  used  to  entrain  air  against  the  resistance  imposed 
by  the  negative  cloud/air  density  gradient. 

However,  lateral  spreading  is  driven  solely  by  the  cloud/air 
buoyancy  deficit.  The  influence  of  ambient  diffusion  and  side 
entrainment  is  ignored.  Here,  our  initial  thought  is  that  the 
assumed  vigorous  unstable  mixing  seems  to  run  counter  to  the 
neglect  of  side  entrainment,  even  though  h/R  is  small.  Cloud 
height  is  determined  by  the  rates  of  lateral  spread  and  cloud-top 
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entrainment.  This  allows  the  cloud  to  slump  before  growing. 
However,  Jensen  (1981)  showed  that  though  side  entrainment  hardly 
affects  the  cloud's  horizontal  spread  rate,  it  may  greatly  reduce 
the  initial  slumping  rate. 

Since  the  authors  believe  that  rapid  mixing  due  to  rupture  of  a 
pressurized  tank  is  the  most  likely  scenario,  the  puff  advection 
velocity  is  set  to  equal  the  ambient  wind  sp'^ed  at  cloud  top.  Data 
from  isothermal  releases  at  Thorney  island  showed  less  early 
advection  than  HEAVYPUFF  due  to  the  initially  stationary  puff  o 
inertia.  The  HEAVYPUFF  equations  are  for  convecti/e  turbulence 
which  implies  relatively  rapid  mixing  So  even  if  the  irsitial  puff 
speed  were  other  than  ambient,  HEAVYPUFF  is  not  really  designed  to 
handle  an  initially  stationary,  quiescently  turbulent  puff,  such  as 
those  released  at  Thorney  Island. 

Even  for  cases  of  rapid  mixing  over  flat  terrain,  we  suspect  that 
cloud  speeds  may  differ  from  ambient,  since  in-cloud  and  ambient 
Obukhov  lengths,  and  L,  will  differ.  I.e.,  cloud/ambient 
differences  in  stability  should  generate  different  vertical 
profiles  of  wind  speed  resulting  in  cloud-top  shear.  The  computed 
speed  difference  might  provide  a  basis  for  a  more  definite  and 
larger  value  for  shear  induced  cloud-top  entrainment  than  the  form, 
Cu»^/h,  now  assumed  in  HEAVYPUFF 's  tke  equation.  As  diffusion 
proceeds  and  cloud-top  shear  and  surface  heat  flux  decay,  use  of 
such  an  algorithm  should  allow  L^,  values  and  hence  cloud-top  speeds 
to  approach  ambient  levels  smoothly. 

Down-slope  flow  is  neglected;  thus,  so  is  lead  edge  entrainment. 
This  probably  makes  the  current  model  unsuitable,  since 
Vandenberg's  Hypergolic  Stockpile  and  Storage  facility  sits  on  a 
mesa  often  just  upwind  of  the  steep  slope  to  the  Santa  Ynez  River 
Valley.  Jensen  (1981)  has  suggested  a  number  of  first  order 
adjustments  for  terrain  slope,  side  entrainment,  building  wake 
effects,  above  ground  releases,  and  the  modeling  of  dense  plumes 
rather  than  puffs,  all  of  which  may  be  suitable  for  eventual 
inclusion  in  HEAVYPUFF.  More  recently,  Kukkonen  and  Nikmo  (1992) 
have  installed  down  slope  components  in  a  cylindrical  box  model. 
The  model's  simplicity  allows  for  fast  runtimes  and  coupling  to 
RIMPUFF  should  also  be  fairly  straight  forward.  However,  source 
term  mechanics  and  chemistry  and  engineering/probability  estimates 
for  flashing/non-flashing  scenarios  are  still  required. 
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(MATHEW/ADPIC  M.  Dickerson,  R.  Lange,  D.L.  Ermak, 
Lawrence  Livermore  Labs) 


MATHEW:  (M.  Dickerson,  C.S.  Sherman) 

MATHEW  is  a  vrind  flow  model  used  to  drive  ADPIC.  MATHEW  makes 
minimal  adjustments  to  an  objectively  analyzed  wind  field  by 
variationally  minimizing  the  3-D  divergence  to  ensure  mass 
consistency.  Unlike  LINCOM,  MATHEW  does  not  include  momentum 
balance  in  its  interpolative  adjustment.  The  description  given  in 
the  MATHEW/ADPIC  user's  guide  is  highly  mathematical  and  optimized 
for  vector  computer  performance.  But  the  main  concepts  can  be 
taken  from  an  earlier  description  by  Sherman  (1978) .  In  essence 
the  operative  error  functional  to  be  minimized  is, 


E(u,v,w)  =  |JJ[a;(u-Uo)  +  a\{v-VQ)  +  aliv^-Wfy) 
^  du  ^  dv  ^ 


(13) 


where  u^,  v^,  w^,  are  the  objectively  analyzed  variables;  X(x,y,z) 
is  the  Lagrange  multiplier;  and  the  ttj  are  Gauss  precision  moduli 
taken  to  be  af  s  1/2  The  CTj  are  deviations  of  the  objective 
from  the  adjusted  field.  This  model  is  unlike  WOCSS  or  MACHV7IND, 
where  vertical  divergence  is  parsed  over  layers  defined  by  a  CSH  to 
produce  a  pseudo-2-D  divergence  equation,  as  in  eqns.  (6-12) . 
Unlike  NUATMOS,  the  vertical  Gauss  modulus,  03,  does  not  try  to 
include  the  effect  of  speed-up  over  hills,  but  the  code  does  have 
stability  dependent  divergence  parsing.  At  any  rate,  minimization 
of  eqn.  (13)  requires  the  solution  of 


u  =  u. 


1  d\ 


2a 


2  ~5x 


V  -  V  +  ^ 

"  ~2  7y 


2  a 


w  =  w  +  ^ 


du  ^  dv 
'dx  1y 


(14) 


subject  to  the  x,  y,  and  z  boundary  conditions, 
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X  S  (u)  =  0 


ny  X  S(v)  =0  (15) 

X  S  (w)  =  0  . 


S(  )  gives  the  first  variation  in  the  bracketed  quantity  and  the  n; 
are  the  positive  outward  unit  normals  on  the  boundaries.  X  is 
obtained  by  differentiating  eqns.  (14a, b,  and  c)  and  substituting 
into  eqn.  (14d)  to  give 


d~X  ^ 
dx^ 


d^X 

dz^ 


(16) 


Eqn.  (16)  is  solved  with  the  boundary  conditions  from  eqn.  (15)  and 
the  adjusted  velocity  field  is  obtained  from  eqn.  (13) . 

Aside  from  being  a  purely  mass  conservative  model,  perhaps  the  main 
limitation  of  MATHEW  is  the  approximation  of  terrain  features  as 
rectangular  blocks  at  the  resolution  of  the  Euler ian  grid.  This 
can  lead  to  excessive  flow  blockage  when  the  terrain  is  steep. 
When  used  with  the  dense  gas  version  of  ADPIC,  the  imposed  stair- 
stepped  approximation  to  terrain  slopes  can  also  produce 
distortions  in  cloud  shape,  since  dense  gas  clouds  may  be  on  the 
order  of  one  grid  cell  high  (see  ADPIC  review  below) .  A  change  to 
terrain  following  coordinates  might  be  useful,  provided  all 
conserved  quantities  can  be  maintained. 

ADPIC:  (R.  Lange,  D.L.  Ermak) 

To  depict  materiel  concentrations  the  dense  gas  version  of  ADPIC 
allows  Lagrangian  particles  to  diffuse  in  Monte  Carlo  fashion 
within  a  3-D  Eulerian  grid  composed  of  rectangular  cells.  This  is 
somewhat  different  from  the  original  neutral  ADPIC  which  used  a 
down-gradient  diffusion  model  for  turbulence  to  solve  the  flux 
conservation  equation. 


dC 

Tt 


^  ^'{CTTp) 


0 


(17) 


Here  C  is  concentration  and  is  a  pseudo-transport  velocity. 
includes  a  mean  advection  velocity,  17^  ,  from  MATHEW,  and  a 
diffusion  velocity,  TT^  ,  both  interpolated  to  each  particle's 
position.  In  dense  ADPIC  D^-  is  not  imposed  as  being  proportional 
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to  the  concentration  gradient.  Instead,  it  arises  from  the 
aggregate  of  random  and  independent  particle  motions.  This  point 
is  not  always  clear  and  both  versions  of  ADPIC  are  rather  different 
from  the  above-reviewed  puff  models.  So  it  may  be  useful  to 
analyze  the  main  equations  in  some  detail.  In  discussing  the 
original  ADPIC,  if  we  use  a  vector  identity,  eqn.  (17)  is 
equivalent  to 


dc 

Tt 


TTp'VC  +  CV^TTp  =  0 


(18) 


Since 


TTp  , 


(19) 


we  have. 


dc 

-Jt 


+  ^  *  CV‘  +  TTj)  =  0  . 


(20) 


MATHEW  has  already  arranged  17^  to  be  non-divergent ,  i.e.. 


(21) 


So  eqn.  (20)  becomes 


dc 

TE 


+  XT^-VC  =  -JTj  ‘VC  -  CV’ZTj  . 


(22) 


Reversing  the  above  vector  identity  now  gives. 


dC 

Jt 


ir^-vc  =  -V-(CU^) 


(23) 


In  dense  ADPIC  the  particles  move  independently  and  randomly.  But 
the  original  neutral  ADPIC  assumes  that  diffusion  velocity  is 
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proportional  to  the  material  concentration  gradient  determined  at 
the  grid  points, 


TF^  =  -K-VC/C  , 


(24) 


where  K  is  a  diffusion  coefficient  appropriate  for  the  ambient 
conditions.  The  material  concentration  gradient  is  equated  with 
the  aggregate  particle  density  gradient.  Each  particle's  diffusion 
velocity  is  obtained  by  interpolating  the  values  at  the  grid  to  the 
individual  particle  positions.  In  this  way  eqn.  (23)  becomes, 


dC 

TE 


+  CT^-VC  =  V*  (2CVC) 


(25) 


Using  non-divergence  and  the  vector  identity  again, 

v-(cu'^)  =  cr-vc  +  cv-u;  =  .  (26) 


This  leaves  the  flux  conservation  equation  in  the  pseudo-velocity 
form. 


dc 

Jt 


V- 


Civz  -  I  VC) 


v.(cr7;) 


0 


(27) 


The  particles  move  within  the  cells  but  advection  velocities:  U,  V, 
W,  and  also  K,  ,  and  C  are  defined  only  at  the  grid  points.  C 

is  obtained  for  each  grid  point  from  a  weighted  sum  of  the  mass 
within  the  eight  surrounding  cells.  This  is  done  simply  to  smooth 
stochastically  induced  unevenness  in  the  distribution. 

The  dense  gas  version  of  ADPIC  adds  an  in-cloud  gravity  flow 
component , 


U  =Xr,  ^TTg-g^iz)  ,  (28) 


to  the  ambient  velocity  field,  U^,  supplied  by  MATHEW.  Here  Ug  is 

the  vertically  averaged  perturbation  of  the  horizontal  wind  due  to 
gravity.  g,(z)  is  a  vertical  profile  function.  The  subscript,  i, 
may  be  horizontal,  h,  or  vertical,  v.  Thus,  the  down-slope  flows 
are  simulated  by  deterministic  advective  components  which  include 
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gravity.  For  the  diffusive  component,  rather  than  the  original 
strict  down-gradient  fashion  of  eqn.  (24)  ,  in  dense  ADPIC  the 
adjusted  velocity  field  is  used  to  displace  the  Lagrangian 
particles  stochastically  in  random  Monte  Carlo  style,  so  that, 

*i  “ 

y,  „  =  y,  ^  VAt  -H  R^(At)  (29) 

z,  ^1  =  z,  +  ^At  +  R^(At)  +  v„At  . 


U,  V,  and  W  are  the  local  velocity  components,  interpolated  from 
the  grid  to  the  individual  particle  positions,  and  R^,  Ry,  and  R^ 
are  independent,  normally  distributed  (Gaussian)  random 
displacement  components.  The  R;  are  specified  as  having  zero  means 
and  the  following  mean  square  (<  >)  properties. 


<R;>  =  2K,At 
^;>  =  2KyAt 

=  2K^(z{t))At  +  vlAt^- 


(30) 


The  horizontal  and  vertical  dif fusivities  for  momentum,  Kj,  are 
derived  from  similarity  theories  and  Vj,  =  dK^/dz  is  a  drift  velocity 
correction.  Drift  correction  is  used  to  counter  the  tendency  for 
particles  to  accumulate  in  low  turbulence  regions.  The  problem  is 
that  the  near  presence  of  a  hard  surface  truncates  the  modeled 
turbulence  scales  and  severely  restricts  particle  motion,  producing 
long  residence  times  near  the  ground.  So  particles  tend  to 
accumulate.  In  ADPIC  this  effect  is  implied  by  approaching  zero 
as  z  approaches  zero.  In  reality,  turbulence  tends  to  be 
positively  skewed  near  the  ground.  So  large  but  intermittent 
burst/sweep  turbulence  events  will  at  times  breach  the  quasi- 
laminar  sub-layer  and  eject  particles  which  would  otherwise  be 
trapped  near  the  surface.  Some  models  include  such  skewness 
explicitly  (see  Baerentsen  and  Berkowicz,  1984) .  However,  this 
generates  a  modeling  paradox,  since  particle  reflection  off  the 
surface  should  reverse  any  skewness  and  in  the  aggregate  make  it 
difficult  to  maintain.  So  drift  velocity  is  a  standard  ad  hoc 
method  used  to  avoid  such  complexities. 

ADPIC  equates  drift  velocity  with  the  vertical  momentum  diffusivity 
gradient,  =  dK^/dz.  This  form  of  drift  velocity  is  non-local  but 
computationally  less  expensive  than  the  original  form  proposed  by 
Legg  and  Raupach  (1982)  (installed  in  a  third  version  of  ADPIC) 
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which  requires  that  the  local  gradient  of  vertical  velocity 
variance  be  computed  instead.  Since  dK^/^z  becomes  large  near  the 
surface,  the  general  effect  is  the  same,  to  give  near-surface 
particles  an  extra  upward  kick. 

The  vertical  diffusivity  is  defined  according  to  surface  layer 
similarity  theory  by 


K 


I 


ku ,  z 


1  + 


(31) 


where  in-cloud  and  ambient  Obukhov  lengths,  and  L^,  differ  by 


+  2gk- 


(P  -  Pa) 


(32) 


and  u,  is  friction  velocity.  p  and  represent  the  vertically 
averaged  in-cloud  and  ambient  densities. 

Simple  bulk  transfer  is  vised  to  account  for  the  ground  flux  of 
heat,  Jg,  into  a  cold,  dense  cloud, 

Jg  =  pCpV^iTg  -  T)  .  (33) 


Rather  than  employing  a  more  involved  surface  energy  balance,  T^, 
the  ground  temperature  has  thus  far  been  equated  with  ambient  air 
temperature,  T^.  T  and  Cp  are  local  layer  averaged  cloud 
temperature  and  heat  capacity,  and  is  a  bulk  heat  transfer 
coeff icient . 

In  inspecting  eqns.  (31-32),  we  note  that  standard  similarity 
theory  defines  L.  in  terms  of  the  friction  velocity  and  kinematic 
heat  flux,  ,  i.e.. 


e 


gkw'e 


(34) 


So  in  eqn.  (32)  for  l'  =  0  (neutral  stability)  we  see  that  the 
fractional  density  surplus  in  the  ADPIC  cloud  is  given  by 
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(35) 


P  -  Pn  ^  wW 
p  ku ,  & 


For  cold  releases  due  to  cryogenic  stores  or  pressure  drops 
iol lowing  rupture  of  a  pressurized  tan)c,  vPW  may  be  of  order  1  ms 

'®K  or  larger,  k  0.4,  ambient  u,  may  be  “  O.l  -  0.2ms  ',  and  9  = 
300 °K.  So  the  equations  seem  to  allow  the  cloud  to  become  unstable 
(Lj  <  0)  at  density  differences  of  "  4%  or  larger,  i.e.,  well  before 
reaching  ambient  density.  If  so,  both  the  form  used  for  and  the 
denominator  in  eqn.  (31)  (given  above  only  in  the  standard 
similarity  form  for  stable  conditions)  must  change  to  reflect 
unstable  conditions.  In  turn  this  should  increase  the  vertical 
diffusivity,  expand  the  cloud,  and  speed  its  approach  to  neutral 
gas  behavior.  In  unstable  ambient  conditions  the  feedback  loop 
appears  positive  because  the  cloud  will  also  tend  toward 
instability  as  the  cloud/air  density  difference  decreases.  If  so, 
accurate  timing  of  stable/unstable  transitions  may  be  significant 
generally  in  modeling  surface-heated  cloud  development. 

The  authors  respond  that  the  above  u,  is  actually  an  in-cloud  term 
which  varies  with  L^,  such  that  cloud  instability  is  unlikely;  in 
any  event  cloud  growth  is  not  affected  much  by  in-cloud  stability. 
On  the  other  hand,  HEAVYPUFF  assumes  that  a  surface-heated  cloud  is 
essentially  unstable.  So  a  convective  velocity  like  w,  can  be  used 
to  help  scale  the  cloud-top  entrainment  rate.  Clearly,  the 
difference  in  physical  assumptions  suggests  more  scrutiny  of  both 
models,  as  well  as  the  timing  and  effect  of  real  cloud  stabilities. 

Above  the  surface  layer,  is  computed  by  the  empirical  form. 


=  CU,ze 


(36) 


Where  C  is  an  empirically  derived  coefficient  and  is  boundary 
layer  height.  Equation  (31)  could  be  replaced  by  more  modern, 
stability  dependent  forms,  such  as  those  of  Sorbjan  (1989). 

In  ADPIC  Ug  arises  from  a  balance  between  gravitational  forces, 
dissipative  cloud-top  entrainment,  and  ground  friction,  after 
assuming  steady  state  gravity  flow  and  hydrostatic  balance.  This 
yields,  for  the  gravity  induced  component  of  the  horizontal  wind. 
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(37) 


=  ~^-[0.5gfp  ~  pjh"]  g{p  -  -  piW,  , 


=  '-■^[0.5g(p  -  pjh^]  -  gift  -  Pa)^^,  “  P^f. 


where  is  the  cloud-top  entrainment  rate,  g  is  gravitational 
acceleration,  h  is  cloud  height  above  the  terrain,  H  is  terrain 
height,  and  Uy  is  the  surface  friction  coefficient,  co^  is  taken  to 
be  proportional  to  K,(h)/h.  £ince  the  densities  are  vertically 
averaged,  this  gravity  component  does  not  seem  to  be  height 
dependent . 

The  gravity  induced  component  of  the  vertical  wind  does  appear  to 
be  height  dependent  and  is  given  by, 


'lx 


dz  g^(z) 


(38) 


Save  for  the  gh(z)  and  gy(z)  adjustments,  eqn.  (38)  looks  like  the 
one  for  which  a  detailed  derivation  was  given  in  the  review  of 
WOCSS  and  MACHWIND  (eqn.  4).  But  here  the  gravity  flow  field,  a 
sub-set  of  the  total  velocity  field,  is  assumed  to  be  non-divergent 
in  itself.  We  question  this  and  also  wonder  whether  the 
hydrostatic  balance  assumed  in  eqn.  (37)  applies  to  such  small 
scales.  The  current  cloud-top  entrainment  term  also  seems  •co 
neglect  the  resistance  provided  by  the  negative  density  gradient 
across  the  cloud/air  interface  and  also  assumes  that  the  cloud  is 
stable.  fs  discus.sed  above,  parametrized  tke  budgets  are  an 
alternative  approach  to  cloud-top  entrainment  which  can  be  used  for 
both  stable  and  unstable  conditions. 

Gravity  flow  should  allow  the  initial  slumping  phase  to  occur  for 
isothermal  cases  and  the  Monte  Carlo  diffusion  should  model  the 
side  entrainment.  However,  since  shear  is  currently  neglected,  the 
extra  lead  edge  entrainment  expected  during  down-slope  flow  does 
not  appear  to  be  modeled  in  dense  ADPIC. 

In  ADPIC,  an  energy  deficit,  Cj,  induced  by  the  cloud/air  enthalpy 
(heat  content  at  constant  pressure)  difference  is  used  to  treat  the 
effect  of  ground  heating.  Starting  from  an  initial  value  of  unity, 
at  each  time  step,  Cj  decreases  slightly  by  the  inverse  geometric 
factor,  1/(1  +  X),  where  X  =  ^tV,,/h  .  Since  the  time  step.  At,  and 
heat  transfer  coefficient,  V^,  are  fixed,  this  rate  varies  only 
with  cloud  height,  li.  Since  dense  clouds  are  typically  much  wider 
than  they  are  tall,  ADPIC  assumes  reasonably  that  cloud  expansion 


61 


due  to  surface  heating  occurs  solely  in  the  vertical  direction. 
Thus,  changes  in  cloud  height  and  temperature  due  to  surface 
heating  are  directly  proportional.  This  implies  that  the  energy 
deficit  reduction  rate  is  largely  regulated  by  the  vertically 
averaged  cloud  temperature,  T.  Conversely,  the  vertically  averaged 
energy  deficit  defines  T  which  in  turn  controls  the  vertically 
averaged  density,  p,  wnich  in  turn  helps  regulate  the  ground  heat 
flux  and  gravity  flow. 

Advection/dif fusion  effects  are  included  by  tagging  each  Laqrangian 
particle  with  an  individual  energy  deficit.  Thus,  Monte  Carlo 
diffusion  transports  the  enthalpy  in  all  directions,  while 
aggregate  displacements  affecting  cloud  height  feed  back  to  the 
energy  deficit  reduction  rate.  Without  deeper  analysis  we  cannot 
say  yet  whether  this  energy  deficit  approach  properly  couples  and 
represents  both  dispersion  and  heating  effects. 

Some  further  algebra,  plus  the  approximation  that  thermal  expansion 
occurs  only  vertically,  leads  to  an  equation  for  the  additional 
particle  displacement  due  to  surface  heating, 


As:.. 


T 

-  1) 
1  +  \ 


(39) 


So  for  each  particle  Az^  is  added  at  each  time  step  to  the  vertical 
displacement  obtained  from  eqn.  (29) .  Also  at  each  time  step,  all 
the  dispersive  computations  are  executed  separately  from  the 
thermodynamical  ones  in  order  to  simplify  the  equations. 

ADPIC  is  a  mature  model  in  the  sense  that  many  features  of 
operational  interest  have  been  added  to  the  basic  scheme.  The 
model  includes  graphical  user  interfaces,  dry  and  wet  deposition, 
plume  rise,  and  explosive  cloud  rise.  This  last  feature  may  be  of 
particular  interest  to  Vandenbcrg  and  should  be  analyzed  in  detail 
in  a  later  study. 

Building  wake  effects  are  not  included  yet  in  ADPIC,  but.  model 
structure  does  not  seem  to  preclude  the  necessary  local  adjustments 
to  the  diffusivity  profiles.  The  rectangular  Eulerian  grid  used  in 
MATHEW  creates  a  stair-step  surface  along  slopes  which  can  lead  to 
significant  flow  artifacts  for  low-lying  dense  gases.  An  ADPIC 
upgrade,  scheduled  for  mid-1993,  will  replace  the  stair  stepping 
with  an  improved  piece-wise  continuous,  linear  interpolation  from 
grid  point  to  grid  point.  However,  MATHEW  will  probably  retain  the 
stair-stepping  beyond  1993.  As  with  the  other  models,  additional 
source  term  chemistry  and  mechanics  are  also  needed  to  determine 
proper  input  values  and  potential  scenarios  such  as  flashing. 
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The  model  compared  quite  closely  with  results  from  two  cases  of 
111.7  “K  (boiling  point)  methane  released  into  an  atmosphere  at  298 
°K.  ADPIC's  equations  have  exact  analytic  solutions  for  these  two 
idealized  cases.  In  both  cases  advection  occurred  only  in  the  x 
direction.  The  first  case  was  isothermal  and  only  vertical 
diffusion  was  allowed.  Vertical  diffusion  was  turned  off  in  the 
second  case,  while  surface  heating  was  turned  on.  This  indicates 
that  numerical  errors  and  errors  due  to  the  finite  number  of 
particles  (5000)  were  quite  minimal  ("1%)  over  the  time  span  of  the 
comparisons  (36  seconds) . 

Further  comparisons  against  experimental  data  would  be  required  to 
analyze  the  quality  of  ADPIC's  new  physics.  In  particular,  we 
might  suggest  that  ADPIC's  thermodynamics  be  tested  under  more 
controlled  conditions  than  can  be  achieved  in  the  atmosphere.  This 
might  be  possible  by  making  comparisons  against  heated  tank  data 
involving  two  miscible  fluids  with  different  densiti.es. 

Execution  time  is  a  remaining  question.  For  a  demonstration  51  x 
51  X  15  grid  over  a  35  km  domain  the  objective  analysis  routine, 
MEDIC,  executes  together  with  MATHEW  in  less  than  one  minute  on  a 
DEC  5000/200  workstation.  In  the  same  demo  five  hours  of  diffusion 
using  K-theory  ADPIC  and  5000  particles  took  less  than  one  minute. 
Monte  Carlo  versions  should  execute  faster,  since  the  aggregate 
particle  concentration  gradient  does  not  need  to  be  e^/aluated  to 
obtain  the  grid  point  diffusion  velocities  which  must  then  be 
interpolated  back  to  the  individual  particle  positions.  Also,  grid 
cells  lacking  particles  can  be  ignored.  So  the  main  computational 
expense  for  Monte  Carlo  ADPIC  is  random  number  generation  which  is 
in  the  process  of  being  made  more  efficient.  This  type  of  random 
displacement  model  appears  to  be  considerably  faster  than  the 
random  velocity  models  used  in  most  Lagrangian  particle 
simulations.  So  the  inclusion  of  ADPIC  within  an  optimal 
diagnostic  modeling  suite  is  tempting,  since  it  is  designed  to 
handle  many  of  the  dense  gas  scenarios  one  expects  for  large 
accidental  releases  at  Vandenberg. 

SUMMARY: 

The  new  ADPIC  is  a  quasi-3-L  model  which  may  be  suitable  for 
emergency  response  use  at  Vandenberg  for  puffs  and  plumes.  It 
employs  random  Monte  Carlo  rather  than  down-gradient  particle 
motions  to  model  diffusion.  Gravitational  terms  are  added  to  treat 
dense  gas  behavior.  ADPIC  is  designed  to  treat  surface  heating, 
cloud-top  entrainment,  and  down-slope  flow. 

For  surface  heated  dense  gases  the  possible  transition  from  stabJe 
to  unstable  cloud  conditions  may  be  of  significance,  since  once  the 
cloud  is  unstable  the  heating  feedback  loop  becomes  positive  and 
mixing  may  increase  considerably.  We  seek  clarification  of  the 
non-divergence  and  hydrostatic  assumptions  used  to  obtain  the 
gravity  flow  terms  and  also  suggest  an  alternative  tke  approach  to 
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modeling  the  cloud-top  entrainment.  Extra  lead-edge  entrainment 
during  down-slope  flow  does  not  seem  to  be  included.  The 
Lagrangian  particles  carry  not  only  material  mass  but  also 
individual  energy  deficits  which  are  used  to  determine  temperature 
changes  and  thermal  expansion  effects.  We  cannot  say  yet  whether 
this  energy  deficit  approach  properly  couples  and  represents  both 
dispersion  and  heating  effects. 

Stair-stepped  surfacing  along  slopes  will  be  retained  in  MATHEW  in 
the  near  term  but  is  to  be  replaced  shortly  in  ADPIC  by  straight 
line  interpolation  between  grid  points.  This  should  be  of 
considerable  benefit  to  low-lying  dense  gas  simulations.  Building 
wake  effects  are  not  included,  but  the  model  structure  does  not 
seem  to  preclude  the  required  adjustments.  Comparisons  with 
analytical  results  indicate  accurate  numerical  execution  for  two 
idealized  cases.  Further  tests  against  relevant  experimental  data 
may  be  required  to  assess  ADPIC 's  underlying  physics.  Source  te  m 
chemistry  and  mechanics  are  also  needed  to  implement  ADPIC 
operationally  at  Vandenberg.  In  general  some  of  dense  ADPIC 's 
physics  may  not  be  fully  mature  or  at  least  remains  to  be  validated 
at  this  juncture.  However,  Monte  Carlo  deployment  of  thousands  of 
Lagrangian  particles  in  a  manner  which  includes  dense  gas  effects 
and  is  fast  enough  for  emergency  response  makes  ADPIC  a  strong 
candidate  for  installation  at  Vandenberg.  Like  most  diffusion 
models  ADPIC  can  probably  be  decoupled  fairly  readily  from  MATHEW. 
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OVERALL  SUMMARY  COMMENTS  ON  THE  MODELS: 


This  summary  is  predicated  on  ready  decoupling  of  the  windflow  and 
diffusion  components  of  the  modeling  systems  reviewed  above.  Based 
on  the  currently  available  documentation  and  perhaps  in  the 
following  order,  dense  ADPIC,  RIMPUFF,  and  CALMET  appear  to  be  the 
most  viable  diffusion  models  for  Vandenberg  applications  at  this 
time.  ADPIC  has  the  large  advantage  that  it  treats  dense  gases  in 
Lagrangian  particle  fashion  and  yet  seems  fast  enough  for  emergency 
response  on  current  workstations.  It  also  contains  an  explosive 
cloud  rise  module.  At  this  point  we  retain  some  reservations 
regarding  ADPIC 's  physics.  Among  the  puff  models  CALPUFF  seems  to 
have  the  most  inclusive  features,  while  RIMPUFF  has  perhaps  a  more 
fundamentally  sound  basis.  Changes  would  be  required  for  either 
puff  model  to  perform  satisfactorily  in  the  Vandenberg  setting. 
Ad Utions  such  as  a  graphical  user  interface  and  output,  spectral 
filtering,  relevant  chemistry,  source  term  specification,  and  dense 
gas  effects  would  be  required  in  CALPUFF.  Chemistry,  source  terms, 
dense  gas,  improved  vertical  diffusion,  and  near-field 
characterization  of  items  such  as  building  down-wash  would  be 
required  for  both  ADPIC  and  RIMPUFF.  Some  features  might  be 
provided  by  modules  already  contained  in  the  ADAM,  /iFTOX, 
HEAVYPUFF,  and  ADPIC  models.  ADPIC  and  RIMPUFF  already  have  UNIX 
based  graphical  user  interfaces.  Coupled  to  RIMPUFF,  HEAVYPUFF,  a 
simple  cylindrical  box  model,  may  also  be  suitable  for  dense  gas 
dispersion  in  relatively  flat  terrain,  where  the  surface  is  heated, 
but  would  need  adjustments  for  slope  flow  which  are  apparently 
manageable.  All  reviewed  diffusion  models  could  profit  by  the 
inclusion  of  diffusion  and  flow  parametrizations  specifically 
addressed  to  changes  across  the  coastline  and  cloud/clear 
interfaces. 

There  is  no  undisputed  winner  among  the  diagnostic  windflow  models. 
A  selection  scheme  using  LINCOM  for  unstable  to  slightly  stable 
conditions  and  MACHWIND  for  stable  conditions  may  be  most  robust, 
as  well  as  cost  effective,  since  these  two  models  are  probably 
nearest  to  operational  form  for  Vandenberg  purposes. 

Taken  individually,  each  flow  model  displays  limitations  which 
leave  cause  for  concern.  Though  in  rough  agreement  with  laboratory 
results,  the  physical  basis  of  the  Froude  number  criterion  and 
models  based  on  them  such  as  WOCSS,  CALMET,  and  MELSAR  are  suspect. 
At  its  present  grid  resolution,  sub-grid  parametrizations  and 
speed,  CALMET  seams  intended  for  regional  rather  than  Vandenberg 
scales.  According  to  preliminary  results  from  CARB's  Los  Angeles 
Air  Basin  study,  CALMET's  windflow  may  deviate  substantially  from 
tower  indicated  winds  when  its  slope  flow  adjustment  scheme  has 
difficulty  defining  a  mean  background  wind  in  complex  terrain. 
CALMET  may  also  not  be  fast  enough  for  emergency  response  on 
workstation  class  computers  at  the  required  500  meter  resolution 
over  a  50x50  kilometer  domain,  while  MELSAR,  LINCOM,  WOCSS,  MATHEW, 
and  MACHWIND  clearly  are.  However,  the  MELSAR  and  MATHEW  windflow 
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schemes  seem  rather  dated.  According  to  available  documentation 
MELSAR's  maximum  resolution  over  a  50  x  50  kilometer  domain  would 
be  1  kilometer,  i.e.,  insufficient  resolution.  MATHEW  also  uses 
problematic  rectangular  rather  than  terrain  following  coordinates. 
The  flow  schemes  in  LINCOM,  WOCSS,  and  MACHWIND  are  essentially 
terrain  following.  Information  concerning  the  speed  of  NUATMOS  has 
not  been  made  available.  We  understand  that  NUATMOS  was  included 
in  the  test  suite  of  models  over  Rocky  Flats  but  documentation  h?s 
not  been  received.  In  general  the  presumption  that  the  background 
mean  flow  can  be  defined  in  the  context  of  complex  terrain  is  more 
viable  over  the  Vandenberg  domain  than  at  larger  regional  scales. 
At  this  time  LINCOM 's  major  drawback  is  that  its  momentum  budget 
treats  only  neutral  flow  physics.  However,  it  is  the  only  windflow 
model  among  those  seriously  reviewed  which  includes  momentum  as 
well  as  mass  conservation  equations. 

Among  the  larger  issues  such  as  availability,  portability, 
maintenance,  validation,  and  total  life-cycle  cost,  LINCOM/ RIMPUFF 
seems  to  be  the  overall  least  expensive  choice,  since  it  has 
already  been  tested  on  a  relevant  Vandenberg  data  set.  User 
support  and  the  availability  of  future  upgrades  is  also  important 
and  contingent  upon  the  stability  of  support  staff.  Regarding 
these  issues,  CALMET/CALPUFF  is  available  in  its  entirety  from  the 
EPA  electronic  bulletin  board.  WOCSS/Adaptive  plume  remains  in  the 
public  domain.  MELSAR/MESOI  is  presumed  available  from  PG&E,  as  is 
NUATMOS /CITPUFF  from  the  EPA.  CALMET/CALPUFF  and  MELSAR/MESOI  are 
supported  by  large  public  agencies.  MATHEW/ADPIC  is  maintained  by 
a  large  staff  at  Lawrence  Livermore  Labs.  WOCSS/Adaptive  plume  is 
supported  by  SRI  International  but  on  a  lesser  scale.  MACHWINU  is 
maintained  by  U.S.  Army  Atmospheric  Sciences  Lab.  LINCOM/RIMPUFF 
has  been  supported  by  the  Naval  Postgraduate  School  and  RISO 
National  Laboratory  of  Denmark.  The  current  status  of  NUATMOS/ 
CITPUFF  has  not  been  ascertained.  Real  installation  costs  and 
short  term  technical  support  are  important  issues  which  cannot  be 
resolved  prior  to  the  fact.  However,  MACHWIND  and  LINCOM/RIMPUFF 
have  both  been  run  on  Vandenberg  terrain  and  conditions. 

It  is  known  that  versions  of  all  the  above  models,  perhaps  save  for 
NUATMOS /CITPUFF,  have  been  run  on  DEC  computer  systems.  So  MARSS 
compatibility  is  not  expected  to  be  a  problem,  provided  hardware 
upgrades  beyond  the  MicroVAX  II  (1  MIPS)  level  become  available  to 
Vandenberg.  At  this  point  it  is  difficult  to  imagine  an  over¬ 
riding  reason  for  operational  use  of  a  good  diagnostic  dispersion 
model  on  less  than  a  25  MIPS  machine,  i.e.,  a  486/50  PC.  So  chis 
review  has  not  been  written  with  the  constraint  of  MicroVAX  II  and 
VMS  operating  system  compatibility  in  mind.  If  such  a  constraint 
were  invoked,  then  PGEMS  and  LINCOM/RIMPUFF  would  likely  be  the 
only  two  remaining  viable  candidates.  Of  course,  conversion  to  a 
GUI  different  from  the  2-D  graphics  based  MARSS  and  retraining  of 
FSC  staff  involves  some  time,  effort,  and  expense.  This  might  be 
saved  if  the  current  hardware  were  replaced  by  DEC  Alpha-chip  based 
VMS  syscems.  However,  Yamada  Science  &  Art  is  scheduled  to  deliver 
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the  prognostic  HOTMAC/RAPTAD  modeling  system  to  Vandenberg  on  a 
UNIX  workstation  in  late  1993.  So  conversion  or  at  least  inclusion 
of  UNIX  based  systems  appears  to  be  inevitable.  The  potential  also 
exists  for  faster,  more  accurate  reality  checks  of  dispex'sion 
results,  if  modern  3-D  point --and-click  GUIs  are  installed. 

A  number  of  other  models  were  also  considered.  The  FLOWSTAR  flow 
model  from  Carruthers  and  Hunt's  firm,  Envirosoft,  Cambridge, 
England,  is  among  the  most  interesting,  since  it  treats  complex 
terrain  physics  for  a  variety  of  stabilities,  albeit  in  a  linear 
manner.  However,  due  to  its  use  of  numerical  solutions,  it  does 
not  yet  seem  fast  enough  for  Vandenberg  emergency  response,  even  on 
workstation  class  computers.  Also  briefly  considered  were  TRAC 
from  Rocky  Flats,  Colorado,  the  EPA  CTDM  and  CTDM+  models,  and  the 
EPA  Urban  Airshed  Model  which  is  apparently  a  precursor  to  CALMET. 
HARM  II  from  Oakridge  Labs  is  reputed  to  feature  dense  gas, 
chemistry,  and  source  term  modules.  However,  the  available 
documentation  was  not  adequate  to  assess  its  merit. 

We  have  yet  to  mention  that  none  of  the  considered  models  contains 
modules  for  such  items  as:  blast  overpressure,  back  calculation  of 
source  strength  from  on-site  receptor  data,  or  overall  risk 
assessment.  Apparently,  ACTA  is  working  on  a  risk  assessment 
model.  The  commercial  model,  SAFER,  from  Dupont  also  contains  some 
of  these  modules.  However,  SAFER  source  code,  documentation 
details,  and  modification  rights  are  apparently  not  available,  even 
though  the  working  algorithms  are  all  based  on  public  domain  EPA 
models.  So  this  model  was  not  seriously  considered. 

We  presume  that  this  evaluation  of  the  technical  and  practical 
merit  of  several  models  will  be  considered  a.long  with  others  in 
making  an  initial  determination  of  detailea  study  for  a  few  of 
them.  We  hope  these  comments  prove  useful  in  furthering  such  an 
effort . 
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